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Motivated by experiments observing self-organization of cold atoms in optical cavities we inves- 
tigate the collective dynamics of the associated nonequilibrium Dicke model. The model displays a 
rich semiclassical phase diagram of long time attractors including distinct superradiant fixed points, 
bistable and multistable coexistence phases and regimes of persistent oscillations. We explore the 
intrinsic timescales for reaching these asymptotic states and discuss the implications for finite dura- 
tion experiments. On the basis of a semiclassical analysis of the effective Dicke model we find that 
sweep measurements over 200ms may be required in order to access the asymptotic regime. We 
briefly comment on the corrections that may arise due to quantum fluctuations and states outside 
of the effective two-level Dicke model description. 
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I. INTRODUCTION 

In recent years there has been rapid progress in con- 
trolling the behavior of ultra cold atoms using a wide va- 
riety of optical techniques. This includes confining atoms 
in optical traps and optical lattice potentials in conjunc- 
tion with tremendous advances in laser cooling. More 
recently it has become possible to study the properties 
of Bose-Einstein condensates (BEC) in ultra high finesse 
optical cavities [l[ . Closely related experiments have also 
been performed on novel hybrid systems combining opti- 
cal fibers on atom chips [3-01 • A central feature of these 
experiments is that one may access the strongly coupled 
regime of cavity Quantum Electrodynamics (QED). In 
this regime a large number of atoms, A'^, exchange pho- 
tons many times on the timescale set by cavity leakage. 
This permits the exploration of coherent matter-light in- 
teractions and the observation of the collective \/N split- 
ting of the resulting eigenstates. It also leads to novel 
forms of collective dynamics and cavity optomcchanics 
Moreover, the light leaving the cavity provides 
valuable information on strongly correlated phases 
[T^ , thereby fostering links between contemporary prob- 
lems in cold atomic gases, quantum optics and condensed 
matter physics. These systems also offer exciting possi- 
bilities as matter-light interfaces for quantum informa- 
tion processing. This wealth of activity is further stim- 
ulated by pioneering circuit QED experiments 1^, 13 . 
which include direct observations of Berry' pha ses ISf . 
vacuum fluctuations [l^ , collective behavior [ItI [isj and 
three-qubit entanglement [l^. 

An important aspect of these developments is the po- 
tential for novel phases and phase transitions induced 
by the cavity light field. The latter mediates long range 
interactions between the atoms which may strongly in- 
fluence their behavior. It was recently predicted that an 
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FIG. 1. (color online). Experimental setup showing cold 
atoms (red) in an optical cavity with transverse pumping [20l - 
[2^. (a) Below the threshold pump power only the pump mode 
is present, (b) Above the threshold the atoms self-organize 
into a checkerboard lattice and are trapped in the interference 
pattern of the pump and cavity beams. The self-organization 
transition for a BEC is described by the onset of superradi- 
ance in an effective nonequilibrium Dicke model. 



atomic cloud with additional transverse pumping under- 
goes a self-organization transition to a spatially modu- 
lated phase 1201; see Fig. [1] [!t] This was confirmed ex- 
perimentally by the Vuletic group using thermal clouds in 
an optical cavity plj . Above a critical pumping strength 
the atoms self-organize to form a checkerboard pattern 
as illustrated in Fig. [TJ This dynamically generated lat- 
tice leads to a strong enhancement of the cavity light 
field due to coherent Bragg scattering. Heterodyne mea- 
surements of the phase of the cavity output also reveal 
the discrete I2 symmetry breaking of the emergent lat- 
tice. More recently, this self-organization phenomenon 
was investigated experimentally using a BEC in an opti- 
cal cavity [H, H^l ■ In this setup spontaneous sublattice 
symmetry breaking coexists with superfiuid phase coher- 
ence giving rise to a novel form of supersolid [23 - [26| . In 
addition it was pointed out that this self-organization 
transition is a dynamical version of the superradiance 
transition in the Dicke model [H [H, H^. The Dicke 
model [28| - l3^ has a long history and describes two-level 
systems or "spins" uniformly coupled to light. When the 
matter-light coupling exceeds a critical value the Dicke 
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model exhibits a continuous phase transition to a state 
with a non- vanishing photon population and discrete par- 
ity symmetry breaking; for a review of the Dicke model 
and its applications in quantum optics see Ref. [ssj . 

In the present cold atoms setting the effective Dicke 
model spin states are two distinct momentum states of 
the BEC Their splitting is therefore controlled 

by the atomic recoil energy, and this enables the Dicke 
model transition to be observed using light with opti- 
cal frequencies. This approach is a close analogue of an 
elegant theoretical proposal by Dimer et al [s^] , for real- 
izing the Dicke model transition using a Raman pumping 
scheme between distinct hyperfine states. These atomic 
experiments also provide a direct implementation of a 
Dicke model Hamiltonian without any additional dia- 
magnetic terms. This circumvents the usual no- go t heo- 
rems for observing the superradiance transition [35|-|37|. 
The experiments also have close connections to work on 
the Collective Atomic Recoil Laser (CARL) [H-liOl. For 
further work on self-organized matter-light systems and 
the possibility of novel phases and phase transitions in 
multimode cavities sec Refs. 14114491. 



A crucial feature of the cavity superradiance experi- 
ments is that they are intrinsically open systems with 
strong pumping and large cavity loss rates [l^, HI] . Any 
account of their physical properties therefore requires a 
nonequilibrium approach. Motivated by this experimen- 
tal situation we recently explored key aspects of the col- 
lective dynamics of BECs in optical cavities [s^l • On the 
basis of a semiclassical analysis of the generalized Dicke 
model presented in Refs. [22|, [23| we obtained a surpris- 
ingly rich phase diagram of nonequilibrium phases and 
phase transitions. Most interestingly, we have found that 
the open system displays two significant features. First, 
for the parameters used in the recent experiments we find 
additional attractors of the long time dynamics that have 
not yet been seen in experiment. In particular, the ex- 
periment suggests a normal state without photons, in a 
region where the semiclassical analysis predicts that the 
normal state is unstable. Second, we find a rich array of 
new phases in experimentally unexplored regions of the 
phase diagram. This includes novel coexistence phases 
and regimes of persistent oscillations. The aim of this 
present manuscript is to shed further light on these per- 
tinent issues and to develop deeper links between theory 
and experiment. To this end we explore the semiclassical 
collective dynamics with a specific emphasis on the emer- 
gent timescales and the observability of the characteristic 
features. A key finding is that these timescales vary quite 
considerably throughout the phase diagram. Under the 
assumption that the effective Dicke model fully describes 
the experimental system [23| our primary conclusion 
is that the O(lOms) duration of the current experiments 
may not be sufficient to reach the long-time asymptote 
in all cases. We discuss the prospects for observing the 
predicted asymptotic states in longer duration experi- 
ments and in other realizations of the nonequilibrium 
Dicke model. We also discuss the possible role of quan- 



tum fluctuations and states outside of the effective Dicke 
model description. 

The layout of this paper is as follows. In Sec jHlwe 
provide an introduction to the recent experiments [22 J23| 
and the associated Dicke model. In Sec. lIIII we discuss the 
semiclassical dynamics of this inherently open nonequi- 
librium system. In Sec. IIVI we present the dynamical 
phase diagram for the presently available experimental 
parameters, and we discuss the nature of the long time 
attractors and the associated time evolution. In Sec. |V] 
we investigate the characteristic timescales governing the 
initial and final stages of the collective dynamics, and we 
discuss the implications for finite duration experiments. 
In Sec. lVII we investigate the phase diagram for a broader 
range of parameters and we discuss the appearance of 
persistent oscillations. In Sec. IVIII we examine the ef- 
fects of contributions which go beyond the effective Dicke 
model and its semiclassical treatment. We conclude in 
Sec. I Villi and provide directions for further research. We 
also include technical appendices addressing the deriva- 
tion of the effective Dicke model, and the location of fixed 
points and their linear stability properties. We also pro- 
vide further details on the phase diagram. In order to 
make the manuscript self-contained we incorporate some 
of the principal findings of our previous Letter [s^] ■ 



II. EXPERIMENT AND NONEQUILIBRIUM 
GENERALIZED DICKE MODEL 



The experiment of Ref. [H consists of a ^^Rb BEC 
with approximately N — 10^ atoms prepared in their 
motional ground state with \kx,kz) = |0, 0). The atoms 
are placed in an ultra high finesse optical cavity of length 
178 /im and cavity loss rate k = 8.1 MHz. As shown in 
Fig. [l]the BEC is subjected to a transverse pump beam 
with Rabi frequency fip, wavevector k and frequency Up. 
The latter is far detuned from the atomic transition fre- 
quency Wa, in order to avoid population of this excited 
level. One may therefore neglect the effects of sponta- 
neous emission. However, the pump frequency is near 
detuned to the cavity frequency Wc, resulting in efficient 
scattering from the pump beam into the cavity and vice 
versa. The coupling strength of a single atom to the cav- 
ity mode is denoted by go and the corresponding level 
scheme is shown in Fig. [5] The experiment is a close 
analogue of a theoretical proposal by Dimer et al [34 1 
for realizing the Dicke model by using Raman pumping 
to couple to two ground state hyperfine levels. A no- 
table difference is that the present experiment exploits 
a Rayleigh scheme, involving distinct momentum states 
rather than internal hyperfine states. This generically 
leads to the presence of a back-reaction term, discussed 
below, which may be avoided in the proposal of Ref. [33 |. 

Absorption and emission of photons yields an effec- 
tive two-level "spin" system (22l . [27| where spin down 
corresponds to the ground state |0, 0), and spin up cor- 
responds to the excited momentum state | ± fc, ±fc) = 



3 



±fe,0)' 



±2k,±ky 




±k,±k) 



FIG. 2. (color online). Level scheme corresponding to the 
experimental setup shown in Fig. [T] The pump beam has 
frequency ujp, wavevector k, and Rabi frequency ilp. The 
strength of the cavity coupling is given by go. The frequen- 
cies of the cavity and the atomic transition are denoted by uJc 
and iUa respectively. Here \kx,kz) axe momentum states of the 
atoms in the BEG, and excited electronic states are denoted 
by a prime. The atomic ground state with \kx, kz) = |0, 0) and 
the symmetric superposition labelled as |±fc, ±fc) constitute an 
effective two-level system governed by an effective nonequilib- 
rium Dicke model. The corresponding level splitting is given 
by Wo = IuJt, where uir = f?]^ jlm is the atomic recoil energy 
resulting from the absorption or emission of a single photon. 
In general, multi-photon processes are required in order to 
couple the ground state to higher momentum states. 



i |afc, /3fc). The latter denotes the symmetric 

superposition of momentum states resulting from two- 
photon emission and absorption processes with k^^k^ e 
{±fc}. In this basis one may introduce collective spin 
raising operators = Y^^=i I fc, ±A;)i i(0, 0| where i 
labels the atoms and is obtained by Hermitian con- 
jugation. The quantum dynamics of this inherently open 
system, with a large cavity loss rate k, can be described 
by the density matrix equation in Lindblad form (Flj 

dtp = ~i[H, p\~ K, [i^^ipp — + ptp^'tp) , (1) 



and Wo = 2uJr, where uJr = h?k^/2'm is the atomic re- 
coil energy. The term involving U = — (l/4)f/g/(W(j — Wc) 
describes the back-reaction of the cavity light field on 
the BEC, and may be interpreted as the AC Stark shift 
due to the appearance of a weak optical lattice in the 
cavity. In the experiment (22j both the pump and the 
cavity are red detuned from the atomic transition, so 
U is negative. However, both signs of U are physi- 
cally achievable. In the atomic ground state, the ef- 
fective cavity frequency Wcg = cu + USz is given by 
Woff = uj - UN/2 = UJc - ujp - Ngl/2{ijJa - ^c) HI, in 
agreement with Ref. [l^. The Hamiltonian © contains 
both CO- and counter-rotating matter-light couplings de- 
noted by g and g' respectively. In the large atom-pump 
detuning limit relevant to the experiment [22| one may 
write g = g' = gi)Vtp/2{ujp - Wq). 

For the experimental parameters used in Rcf. p^ . 
Wo = 0.047MHz and UN = -6.5k/4 = -13.3MHz, 
where the latter is inferred from the observed disper- 
sive shift of the cavity frequency, Weff = uJc ~ uip + 2U N 
[5^ . In the subsequent discussion we will endevor to 
place these experiments in a broader context, and with- 
out loss of generality we will approximate these condi- 
tions as Wo « 0.05MHz and UN « -lOMHz; note that 
the latter differs from the value taken in our previous 
Letter [s^l due to a small discrepancy in the reported 
Hamiltonian in Rcf. (2^. Specifically, the Hamiltonian 
given in Eq. © differs from Eq. (4) of Ref. [HI due to 
a discrepancy in the indicated matrix element M = 3/4 
in the notation of Ref. [22|. This does not affect the 
location of the reported superradiance transition, but is 
important for establishing the broader phase diagram. In 
addition to the energy and timcscalcs appearing in the 
model described by Eqs. ([T]) and ([2]), there is a limit on 
the duration of current experiments which is set by the 
rate of atom loss. In the initial experiments [22| this was 
of the order of 100ms, but is notably longer in subsequent 
experiments j23| . 



where p is the system density matrix, ■0 is the cavity 
photon mode annihilation operator, and t denotes time 
[52 |. The effective Hamiltonian, H takes the form of a 
generalized Dicke model with [H IH [111 



H = ui^U' + ^o^z + us.^jH' + gii^^s- + ^S+) 



(2) 



where S = {Sx,Sy,Sz) is the effective collective spin 
of length N/2 and ^ ± iSy. The derivation of 
Eq. ([2]) from the microscopic description, together with 
a discussion of higher order contributions, is given in Ap- 
pendix |X1 For weak pumping, and in the limit that the 
atom-cavity detuning is much larger than both the pump- 
cavity detuning and the recoil energy [l^, [2^ , the coef- 
ficients are given by w = Wc — Wp — N{5/8)gQ/{uJa ~ ^c) 



III. SEMICLASSICAL DYNAMICS OF THE 
OPEN SYSTEM 



Having discussed the effective Hamiltonian and the 
density matrix equation of motion in Sec. [Hi we now turn 
to discuss the dynamics arising from this model. This is 
essential in order to interpret time dependent nonequilib- 
rium experiments performed in an open cavity (22l . [23| . In 
view of the large number of atoms comprising the Dicke 
model spin states, we will first consider the semiclassical 
limit of this dynamics [s^]. In Sec. IVIII we will briefly 
comment on the role of quantum fluctuations. 



A. Equations of Motion and Symmetries 

The semiclassical equations of motion for the open sys- 
tem described by Eqs. ([1]) and ([2]) are given by 

S- = -tii^o + U\4f)S- + 2z(.gV> + g'r)S., 
Sz = —igt/jS^ + igt/j*S^ + ig'tpS^ — ig'tp*S'^ , (3) 
■0 = — [k + i{u! + USz)] ijj — igS^ — ig'S'^ , 

where = Sx i iSy, k is the cavity loss rate, and we 
neglect the effects of atom loss [22|. The Hamiltonian in 
Eq. ([2]) conserves the total length of the collective spin 
since [8^,5*0;] = for a = x,y,z. Likewise, Eq. ^ sat- 
isfies dtS^ = for all k. As such, the dynamics can be 
explored on the Bloch sphere with |S| = N/2. In addi- 
tion to this conservation law there are further discrete 
symmetries. In particular, the equations of motion in 
Eq. ([3]) are invariant under the parity transformation 

V-^-^, (4) 

as in the equilibrium Dicke model. This symmetry is 
spontaneously broken on passing from the normal phase 
with ■)/' = to the superradiant phase with V' 7^ (2^ . 
[2^ . The equations of motion are also invariant under the 
combined variable and parameter change 

S ^ -S, VoV'*, Lo ^ —uj, g g' ■ (5) 

As we shall see in Sec. lIIIBl both the symmetry in Eq. ^ 
and the duality relation in Eq. ([5|) will have a direct mani- 
festation on the Bloch sphere portraits; the attractors are 
related by these discrete transformations. 

In order to get a full understanding of the behavior of 
Eq. ([3]), it is necessary to address two questions. The 
first regards the nature of the long time attractors. The 
second concerns the full time evolution and its connec- 
tion to the asymptotics. The former is important because 
there are fundamental differences between the dynamical 
phase diagram of the open system and the equilibrium 
phase diagram of the Hamiltonian; even for k — > these 
are generically distinct. The latter question is crucial be- 
cause the open cavity experiments have a finite duration 
and may not always reach the long time asymptote. 

For most values of the parameters, the long-time 
asymptotes are steady states and may be identified as 
stable fixed points. That is to say, for these parame- 
ters there are values of S and -0 for which S = and 
■0 = 0, and these steady states are the eventual fate of 
the semiclassical dynamics for all initial conditions. Wc 
will therefore discuss what these steady states are, and 
where possible, give analytical formulae for them. In 
the following sections we will then use this information 
to present dynamical phase diagrams, describing which 
stable fixed points exist for different values of the param- 
eters. We will then address the dynamical bifurcations 
that correspond to the phase boundaries, as well as ad- 
dressing those cases where the long time asymptotes are 
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FIG. 3. (color online). Schematic illustration of the different 
types of behavior displayed by the semiclassical equations of 
motion in Eq. ([3]), for trajectories on the Bloch sphere with 
]Sj = N/2, (a) Evolution from an unstable fixed point (open 
circle) to a stable attractor (closed circle) of the long time dy- 
namics for all initial conditions. Both the stable and unstable 
points have S = and ^ = 0. However, for the stable attrac- 
tor small perturbations decay, while for the unstable fixed 
point fluctuations grow, (b) As for (a) but with a hyperbolic 
flxed point (cross) having one stable and one unstable eigen- 
mode. Paths first approach and then leave the vicinity of this 
hyperbolic point, before eventually reaching a stable attrac- 
tor with S = and ip = 0. (c) Dynamics exhibiting a stable 
limit cycle with S 7^ and 7^ for all initial conditions. 



more complicated than steady states, such as persistent 
oscillations. We illustrate these possibilities in Fig. [3l 

Before discussing the semiclassical equations of motion 
in Eq. ((3)) , let us comment on some known limiting cases 
that have been studied in the literature. With k = 0, 
g' = Q and C/ = 0, these equations correspond to the 
semiclassical dynamics of the equilibrium Dicke model 
without the loss of photons and without counter-rotating 
terms [s^. This model arises in various contexts and 
has recently been discussed in relation to nonequilibrium 
Cooper pairing (56l - [59| . In this setting, the constituent 
fermions are modelled by Anderson pseudo spins [60j . 
and the cavity light field corresponds to the closed 
molecular channel. A key finding of these studies is the 
presence of collective oscillations. This may also be seen 
by exploiting the integrability of the closely related BCS 
(Bardeen-Cooper-Schrieffer) Hamiltonian [s^ [s^ . The 
same equations of motion also apply to polariton con- 
densation and the synchronisation of oscillators [6l|, . 
Quantum corrections to this collective dynamics have 
also been explored in Refs. (gsI. [6^. 

In contrast to the case when (7' = 0, when k = 0, 
g = g' and U = the equilibrium Dicke model is no 
longer intcgrable. Nonetheless, the model is tractable in 
the thermodynamic limit and displays a mean field super- 
radiance transition. Strikingly, the energy levels reveal 
a crossover from Poisson statistics to Wigner statistics 
in the vicinity of the critical coupling (sil . [3^ . This indi- 
cates the onset of quantum chaos, and is accompanied by 
chaotic attractors in the analogous classical dynamics. 

More recently, Dimer et al [l^l have proposed a novel 
scheme for realizing the nonequilibrium Dicke model de- 
scribed by Eqs. ([T]) and ^ with k 7^ 0. The parameters 
in this effective model are readily adjustable and they 
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focus on the particular case with g = g' and U = 0. A 
notable observation is that cavity losses lead to a shift 
of the mean field superradiance transition, in agreement 
with recent experiments with U ^ [22|, [2J, [131 • 

It is evident from this survey of limiting cases that rich 
collective dynamics is expected to emerge for the more 
general system of equations given by Eq. (jS]) , and in open 
cavity experiments. We shed light on this below. 



connected in the broader parameter space with g ^ g' ■ 
Nonetheless, it is important to distinguish between these 
distinct solutions of the steady state equations of motion 
when g = g' ■ We will discuss the experimental conse- 
quences of this distinction in Sees. llVlandlVl 



1. Superradiant A (SRA) Steady States 



B. Fixed Point Attractors 

In order to get a handle on the possible long time 
steady states, we first enumerate the solutions of the 
equations of motion with S = and ■0 = 0. These fixed 
point solutions may either be stable or unstable, and we 
postpone a discussion of their stability properties until 
Sec. nil CI It is readily verified that the normal state (JJ.) 
is always a possible steady state solution with all the 
spins pointing down, Sz = —N/2 and no photons, "0 = 0. 
Likewise, so is the inverted state (fl") with all the spins 
pointing up, Sz = N/2 and no photons, '0 = 0. More 
generally one may look for non-trivial solutions with a 
non- vanishing photon population and a non-trivial mag- 
netization, Sz- To find these configurations we first note 
that a steady state solution satisfying the first equation 
in Eq. ([3]) automatically satisfies the second equation. 
As such Eq. ([3]) reduces to a pair of complex equations. 
Denoting = -01 + *'02 and S^ ^ Sx i iSy one obtains 



[ujo + u{'^l + ii)]{Sx-iSy) - 



2[{g + g')^Pi+iig-g')^P2]Sz 



(6) 



and 



[ti+iiLj+USz)]{iJi+iiJ2) = -iig+g')Sx-ig-g')Sy. (7) 

In general these equations may be difficult to solve ana- 
lytically. However simplifications occur when U = or 
when g = g' ■ We focus here on the latter since the ex- 
periments of Ref. [2^ correspond to g = g' and negative 
U. We discuss the behavior for g 7^ g' in Appendix IB] 
With g = g' the fixed point Eqs. ^ and (O read 



iLUo + U\^P\^)Sx = 2g{i; + i;*)Sz 
(wo + U\4>\^)Sy = 0, 
{lo + USz — iK)ip = —2gSx- 



(8b) 
(8c) 



It follows from Eq. (|8b|) that there are two classes of solu- 
tions depending on whether ^j, = or a;o-t-t/|0p = 0. We 
consider these two classes in turn and refer to the non- 
trivial steady state solutions as superradiant A (SRA) 
and B (SRB) respectively. Assuming ujq > 0, for C/ > 
only the first type of SRA solution may be present. This 
solution corresponds to the familiar superradiant phase 
in the usual Dickc model where U = 0. For U < the 
second type of SRB solution may exist [s^]- As we shall 
discuss below, in general these solutions are continuously 



Equation (jScj) may be rearranged as an equation for 
and substituted into Eq. ([5a|) : 



uJo{u} + USzf+u;oK^ + 4:g^USl = -8g'^iuj + USz)Sz, (9) 

where we have cancelled a factor of 5^; 7^ from both 
sides of the equation. Using the fixed length spin con- 
straint, Sx = N'^/A — Sz one obtains a quadratic equation 
for Sz- This may be solved to yield [501 



S^ 

N 



UN 



± 



/g27V[4cj2 _ (c/iv)2] -uJqUNk^ 



{UNY{ujoUN- 



(10) 



where the accompanying steady state photon population 
follows from Eq. ([5c|) . In general only one of these roots 
corresponds to a physical solution with |S| < N /2- How- 
ever, as we shall discuss in Sec. IIVBI and Appendix |d1 
there are regions of parameter space where both roots of 
Eq. pO)) are supported; see the regions denoted 2SRA in 
Figs. [51 and [171 In addition, there are two possible signs 
for Sx = ±y/N'^/4: — Sz, where the associated sign of 
is determined by Eq. (|5c[) . This sign choice corresponds 
to the parity symmetry in Eq. ((4|) which is spontaneously 
broken at the superradiance transition. 

The critical coupling strength corresponding to the on- 
set of superradiance starting from the normal (JJ.) or in- 
verted state (tf) is obtained by setting S = (0, 0,^A^/2) 
in Eq. ((9|). One obtains 
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aJo(tL)?: + K^) 



(11) 



where uiz^ = =p and aj„ = UN/2. It is readily seen 
from the Hamiltonian ^ that Wip plays the role of an 
effective cavity frequency for the normal and inverted 
states respectively. For the special case where U = Q this 
agrees with the results of Dimer et al [sj] . In the addi- 
tional limit K = 0, Eq. (|lip reproduces the location of the 
superradiancc transition, gy/N = ^JZmU^/2 for the equi- 
librium Dickc model with counter-rotating terms. More 
generally, Eq. ([TT|) gives the onset of the SRA phase in the 
open cavity system with transverse pumping and g = g' -, 
as recently confirmed experimentally |22| . The explicit 
dependence on k of the phase boundary in Eq. ([TT|) em- 
phasizes the open character of the experimental system. 



2- Superradiant B (SRB) Steady States 



For negative U it is evident from Eq. (|8b[) that another 
class of solutions may be obtained if ujo + U\ip\'^ = 0. 
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Equation (|8al) may thus be fulfilled by taking ip to be 
purely imaginary. It then follows from Eq. (|8cp that uj + 
US, = 0. This yields ^ 

where the magnitude of Sy follows from the normaliza- 
tion condition = 7V^/4. In order to obtain real solu- 
tions for 5*^, we require 3^ + < N^/A. This is equiva- 
lent to the condition g > gi, where 

In order to yield \Sz\ < N/2 we require < In 
contrast to the SRA solution which may exist for either 
sign of U depending on the parameters, the functional 
dependence in Eq. (fT2|) clearly indicates that the SRB 
solution only exists for U < 0. In the special case where 
g = gb and Sy = 0, the SRA and SRB solutions coincide. 

In conjunction with both possible signs for S";; = ± | ^j, | , 
Eq. defines four distinct steady states. These divide 
into two pairs of solutions, where the pairs of solutions 
arc related by the discrete parity symmetry in Eq. Q. 
As we shall see in Sec. IIVBI two of these four solutions 
correspond to stable attractors of the long time dynam- 
ics whilst the other two solutions correspond to unstable 
fixed points; see Fig. [6ljd). 

C. Linear Stability of Fixed Points and More 
Exotic Attractors 

In Sec. nil Bl we discussed the possible fixed points of 
the equations of motion with S = and ?/' = 0. Here we 
turn our attention to the linear stability of these fixed 
points as potential candidates for the long time attrac- 
tors. The calculations are most transparent if we consider 
the instability of the normal (-IJ.) and inverted states (tf) 
where = and Sz = respectively. For arbi- 

trary fixed points the approach is readily generalized but 
is algebraically more involved; the details are outlined in 
Appendix O Writing i/' = "00 + '^V' and S~ = + 5S~ 
where i/jq = Sq ~ 0, 5*^ = and substituting 

into Eq. ([3]) one obtains the linearized equations 

5S- = -iuoSS- T igNS^ T ig'NS^*, 
6tl) = —{k + iu!zp)Sip — igSS^ — ig'SS'^ , 

where Sz = 0, Wip = a; =F a;„ and w„ = UN/2. Parame- 
terizing (5^- = ae-"** -h6*e''''* and SS~ = ce-"** -h d*e"'** 
and equating coefficients with the same time dependence 
one obtains algebraic equations for a, &, c and d. The 
self-consistency equations characterize the possible insta- 
bilities and are given by Eqs. (|C1I) and (jC2p . In the case 
where g = g' the frequencies ?y satisfy 

(77^-iu^-K^)(77^-iUo)=F4g^A^iU;piuo-2iK?7(?7^-cjo) = 0. 

(15) 



The dividing line between exponentially growing and de- 
caying fluctuations corresponds to Eq. ([T5|) having real 
solutions for rj. In this case the imaginary part of Eq. (fT5|) 
vanishes when either rj = ot r]^ = ujq. Demanding that 
the real part of Eq. (fTS)) vanishes when 77 = yields 
Eq. (fTT]) . That is to say, the normal and inverted states 
become unstable at precisely the same point as the SRA 
state becomes possible. For g > g^, Eq. ((T5|) has one 
unstable root and in the language of dynamical phase 
transitions this corresponds to a pitchfork bifurcation. 
In the case of frequencies satisfying r]^ — ujq the real part 
of Eq. ([T5|) vanishes when u!^ = 0. This implies that the 
normal and inverted states also become unstable when 
UJ = ±UN/2 respectively. For values of uj beyond these 
points Eq. (fTSj) develops two unstable roots and the dy- 
namical phase transition corresponds to a Hopf bifurca- 
tion. As we shall see in Sees. IIVI and IVIl all of these 
instabilities describe boundaries in the emergent dynam- 
ical phase diagrams shown in Figs. l4l and 1131 

In the above analysis we have outlined the existence 
of various fixed points and briefly discussed their linear 
stability properties. These considerations are essential 
because more than one of these fixed points may exist 
at a given point in parameter space. For example, the 
normal {]}.) and inverted (ft) fixed points always exist, 
possibly as unstable fixed points, even in the presence of 
the superradiant solutions. As we shall see in Sec. IIVI 
there are in fact cases where more than one stable fixed 
point exists at a given point in parameter space. In addi- 
tion to these coexistence phases, where the final state de- 
pends on the initial conditions, it is also possible to find 
regions of parameter space where no stable fixed point 
exists. In these cases the system may be attracted to 
time dependent solutions such as limit cycles, as found 
in other nonlinear dynamical systems. In the remainder 
of this manuscript we search for the complete set of sta- 
ble attractors of the long time dynamics including fixed 
points, bistable and multistable coexistence phases, and 
time dependent trajectories. 



IV. DYNAMICAL PHASE DIAGRAM OF LONG 
TIME ATTRACTORS FOR g ^ g' AND U <0 

In the previous section we gave a brief overview of the 
simplest fixed point attractors and their linear stability 
properties. In this section we build upon these results 
and establish the dynamical phase diagram correspond- 
ing to the semiclassical dynamics in Eq. That is to 
say, we identify which stable long time attractors exist 
at a given point in parameter space. In order to make 
contact with experiment (22l . [23j in this section we re- 
strict our attention to g = g' and U < 0. We begin in 
Sec. lIV Al bv exploring the phase diagram as a function of 
the remaining parameters, g, uj and U, where the value of 
ujo = 0.05MHz is motivated by Ref. In Sec. HVBl we 
then focus on different points in the dynamical phase di- 
agram in order to clearly expose the nature of the under- 
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lying attractors, including their stability properties and 
their locations on the Bloch sphere. In Sec. II V CI we then 
discuss the characteristic time evolution towards the sta- 
ble asymptotic states, and we rationalize these findings 
using linear stability analysis. In some regimes of param- 
eter space, the time evolution can be rather slow, and we 
characterize where this occurs. We discuss the significant 
implications of these regions of long-lived transients for 
finite duration experiments in Sec. [V] 



A. Phase Diagram of Asymptotic Stable Attractors 

As discussed in Sec. |ll] the experiments of Ref. [l^ 
are performed with wq ~ O.OSMHz, k = 8.1MHz, UN ~ 
— lOMHz and g = g' . We therefore summarise the dy- 
namical phase diagram with these parameters. In order 
to provide some orientation, we illustrate how this phase 
diagram relates to that for the open Dicke model with 
[/ = 0, as a function of uj and g. We consider fixed val- 
ues of the feedback term U , starting from the simplest 
possible case with [/ = 0, and decrease this parameter 
through the experimental value; sec Fig. 

In this nonequilibrium setting, each phase is labelled 
according to the complete set of stable long time attrac- 
tors of the semiclassical dynamics given in Eq. ([3|) . That 
is to say, the phase diagram corresponds to starting the 
system in a wide variety of initial conditions and exam- 
ining the totality of stable end points. In this respect, 
the phase boundaries should be thought of as dynam- 
ical phase transitions, which separate distinct regimes 
of asymptotic behavior. In particular, the blue and red 
boundaries in Fig. 2] correspond to the instability of the 
normal (JJ.) and inverted (ff) states respectively, and are 
given by Eq. (jlip . whilst the accompanying horizontal 
segments correspond to Wqr =0. The gold phase bound- 
ary indicates the critical coupling for the onset of SRB 
and is given by Eq. (|13p . It is important to empha- 
size that whilst all of these dynamical phase boundaries 
may be investigated experimentally, not all of them will 
emerge in a given experiment; the relevant phase bound- 
aries are determined by the initial conditions. In partic- 
ular, the experiments performed so far all begin in the 
normal state (JJ.) with no photons [22|, [2^. However, this 
is not a fundamental experimental restriction, and it is 
essential to survey the totality of attractors for all initial 
conditions, before considering particular initial states. 

For U = Q and w > 0, the structure of Fig. S] mirrors 
the equilibrium phase diagram of the Dicke model, hav- 
ing a transition from a phase at low g where only the 
normal state (JJ.) is possible, to a phase where only the 
SRA state occurs. In the terminology of dynamical sys- 
tems this particular dynamical phase transition occurs 
via a pitchfork bifurcation at gy/N = ■\/wo(w^ -I- k'^)/4:Lo 
[s^ ]: a pair of superradiant fixed points emerge when 
the normal state loses stability. This parallels the sit- 
uation in the equilibrium Dicke model where a pair of 
parity related superradiant solutions emerge at a contin- 
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FIG. 4. (color online). Dynamical phase diagram of the stable 
attractors as a function of the cavity frequency ixi and the cou- 
pling g = g' with parameters k = 8.1MHz and ujq = O.OSMHz 
taken from Ref. [22] . The panels represent different val- 
ues of the feedback term U , going from U = Q (top) to 
UN = — 40MHz (bottom). The second panel corresponds to 
the experimental parameters used in Ref. [S^l . Points (a) - (f ) 
marked in the bottom panel correspond to the fixed point il- 
lustrations shown in Fig. [S] The characteristic time evolution 
at points (b), (f) and (g) is given in Figs. [7| and (8] 



uous phase transition. It is notable that as cj — > 0, the 
critical value of g required for superradiance tends to in- 
finity. This is because for g = g' only the real part of 
-0 drives the polarization of the two-level system via the 
collective couphng g{ip + i/;'^){S'^ + S~); a.s uj ^ 0, ip 
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0) (MHz) 0) (MHz) 

FIG. 5. (color online). Vertical slice through the second panel 
of Fig. |4]with UN = — lOMHz, corresponding to the experi- 
mental parameters used in Ref. [2^ . Variation of (a) IV'P) (c) 
ATg{tp) and (e) Sz on passing from the SRA to SRB phases. 
In the vicinity of both of these transitions there is a narrow 
region of bistability denoted as 2SRA, where both of the SRA 
solutions given in Eq. (|10p coexist; see panels (b), (d) and (f) 
for magnified images of the highlighted regions. 



in Fig. m For example, as a result of the effeetive fre- 
quency shifts induced by negative U, there is a region at 
low g where both the normal and inverted states coexist. 
More strikingly, for UN < —2k in this overlap region, 
there is an extension of the region of superradiant phases 
to lower g, so that the SRB fixed point can coexist with 
both the normal and the inverted states; see for example 
the point (d) in Fig. 21 In such a region, there are multi- 
ple possible stable attractors and the ultimate fate of the 
system depends on the initial conditions. In particular, 
this will lead to multistability of the cavity output field. 
In addition, hysteresis can occur in these multistable re- 
gions. For example, increasing to a large value of g, and 
then slowly reducing the value would lead to superradi- 
ant behavior at the point (d) shown in Fig.lU in contrast, 
slowly increasing g from zero would allow the system to 
remain in the normal state for the same final parameters. 

It is evident from the above discussion that the behav- 
ior of the open system is extremely rich and is fundamen- 
tally distinct from the equilibrium case with k = 0. As 
emphasized above and in Ref. [s^, the behavior of the 
open system is controlled by the stable attractors, which 
do not necessarily coincide with the points of minimal 
free energy. As such, there is a crucial distinction be- 
tween the K — >■ limit of the dynamical system, and the 
equilibrium behavior at k = [66j . In order to address 
experiments with an open cavity one must consider time 
dependent dynamics and n ^ 0. In this setting, equilib- 
rium concepts should only be applied with caution. 



becomes purely imaginary as may be seen from Eq. 
In addition, for a; < 0, the open dynamical system shows 
behavior that could not occur in thermal equilibrium; 
the normal state (JJ.) becomes unstable and the inverted 
state (tf) with Sz = N/2 and -0 = is stable instead. 
It is evident from the Hamiltonian in Eq. ([2]) that the 
inverted state is of higher energy than the normal state. 
However, in contrast to the suggestions of Ref. [g^, the 
relevant question for the experimental realization is dy- 
namical stability, as opposed to minimum free energy. 
This behavior is directly confirmed by the duality of the 
equations of motion given in Eq. ([S]), and is refiectcd in 
the uj — > —to inversion symmetry of Fig. |31 

For negative U, the phase boundaries between the nor- 
mal and inverted states and SRA shift to lower and higher 
frequencies respectively, in accordance with Eq. (|lip ; see 
Fig. m This can be interpreted in terms of a state- 
dependent shift of the cavity freque ncy , uj^ = w =F a;„, 
as suggested by the Hamiltonian © [Ifl S^I ■ Within the 
region of overlap the SRB phase may be stabilized, as 
shown in Fig. |4l In particular, this results in a change in 
both the intensity and the phase of the cavity light field 
as indicated in Fig. [5] In addition, in the vicinity of these 
transitions between SRA and SRB, a narrow coexistence 
region emerges, denoted 2SRA, where both solutions of 
Eq. ([TU)) are physical; see Appendix ID] Indeed, such coex- 
istence phases are abundant in the phase diagram shown 



B. Nature of Attractors 

In the previous section we have considered the dynam- 
ical phase diagram with g = g' and U < 0. Here, we 
provide a detailed discussion of the nature of the long 
time attractors indicated in Fig. |31 We focus on the rep- 
resentative points (a) - (f) shown in Fig. |4l and chart 
the associated motion of the spins on the Bloch sphere 
with |S| = N/2. In Fig. [BJ we show the results obtained 
by numerical integration of the differential equations in 
Eq. (j3]) , using an adaptive time step fourth order Runge- 
Kutta routine from the NAG library [g^]. In addition 
to the characteristic trajectories, in Fig. IH] we indicate 
the locations and nature of the various fixed points with 
■0 = and S = 0. The nature of these attractors is deter- 
mined by the number of unstable eigenmodes for small 
fluctuations. If no eigenmodes are unstable, the fixed 
point is stable and is indicated by a filled circle. If one 
eigenmode is unstable it is a hyperbolic fixed point (or 
equivalently a saddle node) and is marked as a cross. If 
two eigenmodes are unstable it is an unstable fixed point 
and is represented by an open circle. For wq <SC t there 
are never more than two unstable eigenmodes, meaning 
that the state of the photon field rapidly comes to follow 
the state of the collective spin. 

In order to gain some orientation we discuss the indi- 
vidual panels in Fig. [51 For the parameters used in panel 
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FIG. 6. (color online). Bloch spheres with |S| — N/2 corre- 
sponding to the points (a) - (f ) in Fig. [4] showing the fixed 
points where S = and "0 = 0. We distinguish between sta- 
ble fixed points (filled circles), unstable fixed points (open 
circles) and hyperbolic points with one stable and one un- 
stable eigenmode (crosses). The presence of more than one 
stable attractor on the Bloch sphere corresponds to a coex- 
istence phase in Fig. [l] The trajectories show the typical 
time evolution, where the initial conditions are chosen in or- 
der to illustrate the attractors. In particular, the timescale 
for approaching the SRA fixed point in panel (b) is signif- 
icantly longer than the timescale for approaching the SRB 
fixed point in (f). In the language of dynamical systems the 
transition from (a) to (b) is a pitchfork bifurcation where a 
single mode goes unstable. The transition from (c) to (a) is 
a subcritical Hopf bifurcation where two modes go unstable 
simultaneously. The transition from (c) to (d) involves the 
appearance of eight additional fixed points; these correspond 
to two stable and two unstable SRB fixed points and four 
hyperbolic SRA fixed points. The transitions from (d) to (e) 
and (e) to (f) are inverse pitchfork bifurcations in which a 
pair of hyperbolic SRA fixed points coalesce at a previously 
stable fixed point. 



(a) there is one stable attractor on the Bloch sphere cor- 
responding to the normal state (-IJ-). This is the ultimate 
fate of the system for all initial conditions as indicated 
in Fig. m the inverted state (ft) corresponds to an unsta- 
ble fixed point. In passing to Fig.|6l^b), the normal state 
(JJ.) becomes an unstable hyperbolic fixed point and an 
SRA attractor with a non-trivial magnetization 5"^ (or its 
parity symmetry partner on the other side of the Bloch 



sphere) governs the long time dynamics. In contrast, 
in panels (c), (d) and (e) we see multiple stable fixed 
points corresponding to the coexistence phases, -IJ- + tTi 
SRB+ 4 tr and SRB+ ]). respectively; see Fig. S] For 
these parameters, the final state of the system depends 
on the initial conditions, and we only highlight some typ- 
ical trajectories. Nonetheless, the totality of stable fixed 
points completely accounts for the possible asymptotic 
behavior and therefore discriminates between different 
dynamical phases. In panel (f) we see both stable and 
unstable non-trivial fixed points corresponding to the su- 
perradiant phase SRB; see Fig. 51 Note that in Fig.[S]we 
have focused on cases with w > 0. For w < the fixed 
points can be immediately found from the duality under 
the transformation given in Eq. ([5|) . This corresponds to 
inverting the Bloch spheres shown in Fig. [5] 

While the attractors determine the long-time asymp- 
totic behavior of the system, there are cases where the 
evolution proceeds on surprisingly long timescales. This 
may be seen in panels (a)-(d) in Fig. |6l where the 
timescale for relaxation towards the fixed point is much 
longer than the period of the orbits encircling it; indeed, 
the latter cannot be resolved in these panels, leading to 
a dense covering of the Bloch sphere. In addition, the 
total time interval for approaching the SRA fixed point 
in Fig. IHb) is much longer than that for approaching the 
SRB fixed point in Fig. [6jf). We will investigate this 
crucial distinction in more detail below. 



C. Time Evolution 

In order to shed light on the dynamical distinction be- 
tween the SRA and SRB fixed points, we examine the 
time evolution starting near the normal state (-||) for the 
parameters used in Figs. [6jb) and (f) with w > 0; see 
Fig. [T] By specifying an initial condition which is not 
a stable fixed point, the dynamical traces correspond to 
a sudden quench into the superradiant state. However, 
since the initial state corresponds to an unstable fixed 
point, the dynamics would remain stuck in the absence 
of noise or quantum fluctuations. In practice, these in- 
herent fluctuations will destabilise the initial state and a 
non-trivial time evolution will take place. 

In order to probe the intrinsic quench dynamics we 
consider two different approaches for perturbing the ini- 
tial condition. The first approach is to displace the initial 
state by Sx = Sy ~ VN, corresponding to the charac- 
teristic size of quantum fiuctuations in the initial state; 
the subsequent semiclassical time evolution gives the pale 
gray trajectories in Fig.[7l The second approach is to use 
a Wigner distributed ensemble of initial conditions in or- 
der to incorporate harmonic fiuctuations around the nor- 
mal state (-U-); see for example Refs. [6l,|6l] and Appendix 
IfI The corresponding time dynamics is represented by 
the black lines in Fig. [71 It is readily seen that both the 
\/N displacement and Wigner approaches are in quan- 
titative agreement regarding the overall timescales for 



10 



0.4 
0.2 


-0.2 
-0.4 

100 

80 
60 
40 
20 




- 03=3(iMHz, gVN=0.9MHz 0.5 

SC 

Wigner _„ g ^ 


A^AAAAAi 

) 1 2- 




100 

50 

, 




M 11 1 1 1 





10 20 30 40 50 60 70 80 
t (ms) 




0.4 - to=10MHz, gVN=0.8MHz 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 
t (ms) 



FIG. 7. Time evolution in the superradiant regime with initial 
conditions that are close to the normal state (JJ-). The top two 
panels (and insets) correspond to point (b) in Figs. |4] and [6l 
The bottom two panels correspond to point (f). In each case, 
the evolution of the semiclassical equations for a single initial 
condition with Sx = Sy = \/N is shown in gray (marked SC). 
The average evolution for a Wigner-distributed ensemble of 
initial conditions is shown in black. The insets in the top two 
panels show magnified images of the highlighted regions. 



evolution towards the SRA and SRB fixed points. Al- 
though the amplitude of the collective oscillations is par- 
tially washed out by the Wigner distribution of initial 
spin states, oscillations of the same period nonetheless 
remain in these examples. 

Comparing the cases shown in Fig. [3 there is clearly 
a significant difference in the relaxation timescales. For 
evolution towards the SRA fixed point shown in the up- 
per two panels of Fig. [71 there are oscillations at frequen- 
cies of a few kHz, with a decay time of order 100ms. 
In contrast, for evolution towards the SRB fixed point 
shown in the bottom two panels of Fig. [71 the timescale 
of the "oscillation" is similar (~ 0.2ms), but only one 



FIG. 8. (color online). Time evolution starting close to the 
normal state (J|) for the parameters at point (g) in Fig. [Jj As 
in Fig. [T] both the semiclassical dynamics for a single start- 
ing point with Sx = Sy = y/N (marked SC), and a Wigner 
distributed set of initial conditions are shown. The latter is 
indicated in black, and the semiclassical trajectory is shaded 
to match the Bloch sphere shown as an inset; for these initial 
conditions the trajectory almost covers the Bloch sphere. 



oscillation occurs before the steady state is reached. As 
discussed in Ref . (soj , the remarkably slow dynamics near 
the SRA fixed point is related to the proximity to a dy- 
namical phase boundary in the extended parameter space 
where g and g' are allowed to vary independently. We will 
return to explore this point further in Sec. IVIBI 

Thus far we have considered the time evolution at 
points (b) and (f) in Fig. [H corresponding to w > 0. 
For completeness we should also consider the time evo- 
lution at point (g) with uj < 0. Owing to the duality 
in Eq. ([5]), the dynamics starting from Sz = ~N/2 with 
w > is related to the dynamics starting from Sz = N/2 
with UJ < 0. We should therefore consider the quench 
dynamics starting from 5*^ — and w < sepa- 

rately. This is illustrated in Fig. [SI along with the ac- 
companying dynamics on the Bloch sphere. It is notable 
that the dynamics in this regime where w_ < has a 
remarkably long timescale. From the Bloch sphere, and 
the time dependence of Sz, it is clear where this comes 
from: the trajectory first spirals around the unstable nor- 
mal state (-11), growing in amplitude, until it reaches the 
stable manifold of the hyperbolic inverted state (ff), from 
which it then transfers to spiral around the stable attrac- 
tor. As such, almost the entire Bloch sphere is covered 
by this trajectory, and a very long waiting time of order 
0.2s is required before reaching the asymptotic state. In 
contrast to Fig. [Tithe long-time asymptote is not reached 
after 100ms. In Sec. [Vlwe will consider the implications 
of this type of behavior for the phase diagram obtained in 
the experiments of Ref. [22| which have a finite duration 
of order 10ms. 
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FIG. 9. (color online). Characteristic timescales of the semi- 
classical dynamics. We use the same parameters as in Fig. |4] 
with UN — — lOMHz corresponding to the experiments of 
Ref. [22|. The top panel shows the time required for the in- 
stability of the initial ground state (|L) to develop; in regions 
where this state is stable this timescale is taken to oo. The 
lower panel shows the asymptotic timescale for approaching 
the final stable state. In the regions of multistability, the rate 
of attraction to the final state for the given initial conditions 
(ll) is shown. In both panels the timescales show a strong 
variation throughout the phase diagram, with notable impli- 
cations for finite duration experiments. A vertical slice along 
the red dashed line is given in Fig. 1101 



V. GROWTH AND DECAY TIMES AND 
IMPLICATIONS FOR FINITE DURATION 
EXPERIMENTS 

In the above discussion we have seen that the timescale 
for reaching the asymptotic attractors varies considerably 
between different points in parameter space. In order to 
make contact with experiment [2^ it is therefore crucial 
to understand how this timescale varies throughout the 
phase diagram. In this section we address this key issue. 

It is evident from Fig. [8] that in order to characterize 
the temporal evolution we require at least two principal 
timescales. The first is the timescale for departing the ini- 
tial state, and the second is the timescale for approaching 
the final asymptotic attractor. Both of these timescales 
may be extracted by linearizing around the initial and fi- 
nal states as appropriate, and calculating the eigenvalues 
(frequencies) using the methods outlined in Sec. lIII Gl and 
Appendix [Cl In Fig. [9] we use these eigenvalues to plot 
the characteristic time for the normal state to become 
unstable, if it does so, and the characteristic decay time 
in the approach towards the final state; in the case of 



asymptotic coexistence phases we focus on the state that 
is actually reached in a quench experiment that starts 
close to the normal state (JJ.). It is clear from Fig. [5] 
that both of these fundamental timescales vary signifi- 
cantly throughout the phase diagram. In particular, in 
the region where uj- < 0, both the timescale for the ini- 
tial dcstabilization of the normal state and the timescale 
for decay towards the asymptotic state are increased. 
The combination of these two timescales provides a lower 
bound on the overall duration of the intrinsic dynamics. 
For a;_ < we therefore expect a slower approach to the 
asymptotic regime. In order to gain a better handle on 
this issue, we provide an analytic discussion of the con- 
stituent timescales below. We begin in Sec. IV AT] with an 
analysis of the initial growth times before examining the 
final asymptotic decay times in Sec. IV A 21 In Sec. IV Bl 
we then consider the implications for experiments which 
monitor the photon intensity over a finite time interval. 

A. Growth and Decay Times 

1. Growth Times 

The initial growth time in Fig. IHJa) can be understood 
and estimated analytically, by using the linearization dis- 
cussed in Sec. IIII CI Considering Eq. for the normal 
state with Sz = —N/2 the eigenvalues ij obey: 

[(?7 + iK)^ - ujl]i^^ -^o)~ 4c^o^-5'iV = 0, (16) 

where aj_ = lu — UN/2. To obtain the growth rate for 
a given set of parameters we must find the solution of 
Eq. p6p with the largest positive imaginary part, 77". 
Solving Eq. (jl6p numerically we plot the corresponding 
growth time l/rj" in Fig.lHUa). It is evident that there are 
distinct growth times in the top and bottom portions of 
Fig.inja) separated by the critical line ui- = 0. The origin 
of this distinction may be traced to a change in behavior 
of the root structure of Eq. (fTB|) . For the parameters 
shown in Fig. E^a), when uj- > a single one of the 
four roots goes unstable, whilst for lu^ < two roots 
with same imaginary parts go unstable simultaneously. 
In the language of dynamical systems, the first scenario 
corresponds to a pitchfork bifurcation, whilst the latter 
corresponds to a Hopf bifurcation [t^I • 

In order to gain a quantitative handle on the observed 
growth times we may exploit the small parameter ujq/k in 
order to find analytic approximations for rj. Anticipating 
that I77I ^ <C K we neglect 77 in the combination ij + in. 
This yields a quadratic equation which provides 

V^±^u:,--^^^±.o^l-^-j , (17) 

where g~ is given in Eq. (jlip and the plus sign corre- 
sponds to the growing mode. For g > g", 77 is imaginary, 
and the characteristic growth rate is indeed of order wq as 
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FIG. 10. (color online). Vertical slice through Fig. [5] with 
gy/N = 0.6MHz. (a) Initial growth time obtained by linear 
stability analysis around the normal (JJ.) state. The open cir- 
cles show the exact semiclassical results obtained numerically 
from the quartic Eq. (|16|) . The solid line corresponds to the 
approximation given in Eq. H18|) where the plus sign corre- 
sponds to the unstable mode, (b) Asymptotic decay time 
obtained by linear stability analysis around the appropriate 
final asymptotic state. The open circles correspond to the 
eigenvalue of [ryl — M| =0 with the largest imaginary part, 
where M is given by Eq. HC5|I . The solid gold (mid gray) 
line corresponds to the imaginary part of Eq. (|20p and the 
solid yellow (light gray) line corresponds to the exact result 
in Eq. (|22p . In both panels, the shaded regions indicate where 
the normal state (JJ.) is stable. Both timescales show a strong 
dependence on the parameters. 



we assumed. However, for w_ < it is evident from the 
first form of Eq. ^7}i that rj E M. and therefore Eq. (fTTj) 
cannot apply in the lower region of Fig. IHJa) . To fully de- 
scribe the observed behavior it is necessary to carry out 
the expansion of rj to higher order in loq/ n. Parameter- 
izing 77 = ±woy^ 1 — [g/gaY + <5±, where 5± <. LOn k., 
and substituting into Eq. (fT6)) one obtains 



?7 « ±ujo\ 1 - ^ 



(Cj2 +^2)2 • 



(18) 



As shown in Fig. Iior a*). the root with the plus sign in 
Eq. (fTSl) accurately reproduces the exact growth times 
observed in Fig.[9l[a). In particular, for ui- > one may 
neglect the last term in Eq. (|18p . but for w_ < it is 
essential to retain this contribution as the leading term 
no longer corresponds to a decay rate. In the former 
case where a;_ > 0, the decay rate is typically of order 
Wo = 0.05MHz; this corresponds to a timescale of 20/is, 



in agreement with the upper region of Fig. |9lja) and the 
central portion of Fig. fTUfa). In the latter case, where 
a;_ < 0, the decay rate is of order ujq/k ~ 0.3kHz, where 
we use the functional dependence of the critical coupling 
in Eq. ([TT|) . This corresponds to a much longer timescale 
of order 3ms, in agreement with the lower region of 
Fig. ^a.) and the left hand side of Fig. rTOT a). In fact 
the exact timescales diverge on approaching the phase 
boundaries where the -U- state is stable; see Fig. [TUTa) . 
This may also be interpreted as critical slowing down [50| . 
For example, approaching the left region of ||. in Fig. fTOl a) 
from the right, one has rj V ga ~ g] this is analogous 
to an equilibrium mean field exponent zv = 1/2. On 
the other hand, approaching the left region of JJ. from the 
left, one has 77 ~ (w — UN/2). This latter behavior arises 
from the second term in Eq. and only exists in the 
open system with k 7^ 0; it is an analogous to a critical 
exponent zv = 1. A divergent growth time may also be 
seen in the lower left corner of Fig. EJa), since 17 — >■ 
in the second term of Eq. (fT5|) . For recent discussions of 
critical behavior in driven open cavities see Refs. (7ll - [7^ . 



2. Decay Times 

Turning to the asymptotic decay time in Fig. (Htb), 
we consider the approach towards the three stable fixed 
points, ft, SRA and SRB which differ from the initial 
state, JJ.. In order to extract the associated timescales we 
must linearize around the asymptotic fixed points and 
find the eigenvalue with the smallest imaginary part. 

For decay towards the inverted state (ft), we may in- 
voke our previous result in Eq. (|18|) with the replacement 
u;_ — > In order to extract the timescale govern- 

ing the approach towards SRA we must linearize around 
this fixed point. Following the general approach used in 
Sec. IIII Cl this yields the characteristic equation 



V 



V 25, 



'^7 



\2gS,-UijS., 



(19) 

where cj + USz and a)o = + t^lV'P S'l'e useful 
variables suggested by Eq. (jSj; for more details see Ap- 
pendix O and Eq. (|C11|) . The ultimate timescale control- 
ling the decay towards the fixed point is governed by the 
slowest roots of Eq. p^ . Anticipating that these have 
\ri\ ^ CjqN /2Sz <Si K we may once again neglect the term 
77 in the combination i] + in. This yields a quadratic 
equation for 77 with solutions r/ = ±770 where 770 € M. 
In order to refine this approximation we parameterize 
•q — ±7]Q + S± where (S± <C <C n, and substitute into 
Eq. (|19p . Retaining terms up to linear order in S one 
obtains 



77 w ±?7o 



2iKU}QUj\2gSz — USxip\'^ 



■:,2 



^2)2^, 



(20) 



In the limit [7 = one finds Im(77) w —KLo^jiuy^ + k^), 
where we use Eq. ^ to substitute for Sz- This agrees 
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with our previous findings jSOj , and yields a characteris- 
tic decay rate of order Wq/k ~ 0.3kHz. This is consistent 
with the 3ms timescale found in Figs.[9l^b) and fTOl b). In 
addition, in the limiting case of the SRA phase where 
■0 = and Sz = —N/2, we recover the characteristic fre- 
quencies of the normal (4) state as given by Eq. (fTS]). As 
shown in Fig. [TUT b). Eq. ([^0)1 accurately reproduces the 
results obtained by direct numerical solution of Eq. (|19p . 

For decay towards SRB, the results are much simpler. 
In this case the characteristic frequencies satisfy 



[?7(?7 + iK) - AgU^^^Syf = 0, 



(21) 



where ip = ipi + iil)2\ see Appendix [Cl In the case of the 
stable SRB fix ed points discussed in Sec. IIIIB 2| where 
Sy) = ±(j-\/ (~<^o/f^), one readily obtains exact 
results for the repeated roots of Eq. (PT|) : 



-I— ± 
2 




(22) 



where we use the fact that 5^ = i/j^n^^g^'^ — g^'^)/4U. 
The decay towards SRB is governed by the slowest mode 
corresponding to the positive root in Eq. (f^ . Since 
Wo <C K we may Taylor expand this root to obtain 
Im(77) 



'2LUQyJ lg\ — 1 + ©(wq/k). As shown in 



Fig. [TUfb). the exact analytic result ^F2^ is in agreement 
with the numerical solution of Eq. pip as required. 

From the above analysis we see that the charactcrstic 
decay rates towards SRA and SRB are of order Wg /k and 
Wo respectively. These correspond to decay times of the 
order of 3ms and 20/xs which yields a faster approach to- 
wards SRB in comparison to SRA. This is confirmed by 
the typical trajectories on the Bloch sphere as shown in 
Figs.El^b) and (f). It is readily seen that many more or- 
bits are executed in reaching the SRA fixed point. We see 
that the SRA and SRB attractors differ in their dynamic 
characteristics, in addition to their steady state forms. 



B. Photon Intensity Map Extracted at 
Intermediate Times and Implications for Finite 
Duration Experiments 

Having examined the characteristic growth and decay 
times across the phase diagram, we now consider the con- 
sequences for finite duration experiments. In particular, 
for w_ < 0, the semiclassical dynamics predicts relatively 
long growth times and comparably long approach times, 
as shown in Fig. [HI In Fig. [TT] we compare the resulting 
photon intensity map obtained after a hypothetical in- 
finite duration experiment, to that obtained after 10ms. 
It is readily seen that the lower region, corresponding to 
w_ < 0, has not fully reached the asymptotic regime. 

In order to make close contact with the experiment of 
Ref. we should also incorporate the details of their 
data aquisition scheme. In particular, the photon in- 
tensity map is obtained by increasing the laser intensity 
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FIG. 11. (color online). Photon intensity maps showing \ip\'^ 
after two different time intervals, with initial conditions that 
are close to the normal state (J|) with Sx ~ Sy = \/N . We 
use the same parameters as in the second panel of Fig. |3] 
with UN = — lOMHz, corresponding to the experiments of 
Ref. [2^. (a) Intensity map obtained in the final asymptotic 
state with t — >■ oo, showing the distinct regions of SRA, SRB 
and 1^. (b) Intensity map obtained after 10ms showing good 
qualitative agreement with the asymptotic attractors in panel 
(a), but with a slower approach in the SRA regions. A cross 
section of panel (a) with g^N = l.OMHz is provided in Fig. [5] 



over a 10ms interval and recording the photon intensity 
during this period. This procedure is then repeated for 
other detunings and an intensity map is generated. In 
order to facilitate a direct comparison, we incorporate 
the effects of the sweep in our numerical simulations, 
where we take g"^ oc t. In Fig. [T2la) we show how for 
one value of w, j'i/'p evolves with increasing matter-light 
coupling where we take g^N = (t/to) x 2.5MHz^, where 
to = 10ms in Figs. HHa) and fTifb), a nd = 200ms in 
Fig- HHc). As discussed in Sec. II V CI we present both a 
single semiclassical trajectory with Sx — Sy = -v/lV, and 
the results of Wigner distributed initial conditions. As 
found previously, it is readily seen from Fig. [T^a) that 
quantum fluctuations reduce the oscillations in the pho- 
ton intensity, but that the overall dependence conforms 
to the semiclassical analysis. Moreover, for this set of pa- 
rameters, the results of the 10ms sweep are also in good 
agreement with the steady state photon intensity. We 
may therefore use the semiclassical approach to map out 
the resultant phase diagram. 

A notable finding is that in other regions of the phase 
diagram, sweeps longer than 10ms may be required in 
order to reveal signatures of the asymptotic attractors. 
To see this more clearly, in Figs. IT^ b) and (c) we show 
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FIG. 12. (color online). 



found by increasing oc t, 



and recording the value achieved as a function of time. We 
use the same parameters as in Fig. [4] with UN = — lOMHz, 
corresponding to the experiments of Ref. [22|. The sweep is 
chosen so that g'^N = (t/to) x 2.5MHz^, where to = 10ms 
in panels (a) and (b) and to = 200ms in panel (c). The top 
two panels correspond to the experimental sweep duration 
used in Ref. [22]. (a) Comparison of the steady state value 
of with that obtained by semiclassical evolution (marked 
SC) and Wigner distributed initial conditions, for the value 
of g\/N reached at a given time with uj — 40MHz. (b) Pho- 
ton intensity map obtained after a 10ms sweep, (c) Photon 
intensity map obtained after a 200ms sweep. In comparing 
these figures to Fig. 5 of Ref. [S^] one should note that the 
vertical scale on our intensity plots is the photon frequency uj 
appearing in Eq. ([2}. In comparison, Ref. [S^l use the detun- 
ing of the pump from the bare-cavity frequency as the vertical 
scale, hence the inverted and shifted axis. 



the final asymptotic state of the semiclassical dynamics 
has not been reached. 



VI. GENERAL PHASE DIAGRAM 

Having discussed the phase diagram and the collective 
dynamics for the experimentally explored case with g = 
g' and U <Q [22| , we now consider the broader parameter 
space. In Sec. IVI Al we consider the case with g = g' and 
[/ > 0, and in Sec. IVI B] we examine the case with g ^ g' . 
A notable feature of both of these cases are parameter 
regimes in which no stable fixed point exists, and for 
which persistent oscillations arise. 



A. Phase diagram for 



9^9 



and [/ > 



The sign of U can be varied by switching between red 
and blue detuning of the cavity light field with respect 
to the atomic transition. This may be seen from the 
derivation of the effective Dicke model Hamiltonian, as 
outlined in Sec. |TT] and Appendix [X] In Fig. [13] we show 
the resulting phase diagram with g = g' and U > 0] this 
is the analog of the phase diagram shown in Fig. [Jj 

With U > 0, the phase boundaries at uj± = shift 
in the opposite direction to those obtained with U < 0. 
As such, the boundaries separate, rather than overlap, as 
shown in Fig. [T31 Instead of finding coexistence phases, 
as found in Fig. [4l a regime of persistent oscillations 
emerges as shown in white in Fig. 1131 In this region, 
no steady state is ever reached, and the photon number 
continues to oscillate periodically at long times. This is 
illustrated in Fig. [T4| 

For these persistent oscillations with g — g' it is possi- 
ble to characterize their behavior analytically. In this 
case the emergent state is a limit cycle [50] ■ To see 
this one may note from Fig. [TJ) that the asymptotic 
behavior has constant Sz, and in fact Sz = —lu/U. 
From the equations of motion in Eq. ^ we see that 
Sz = —ig{'4' + ~ S~). However it also clear from 

Fig. [H that 5y ^ and so S+ ^ . We therefore re- 
quire Re('0) = 0, as found in the SRB steady state. With 
these conditions on Sz and V', Eq. ^ simplifies. Writing 
S^ = re^*^, where = iV^/4 — uP' jU"^ is a constant of 
the motion, one obtains 



the results of the semiclassical evolution in Eq. ([3|) , with 
a/ZV displacement of the initial state, and a sweep profile. 
It is readily seen that for aj_ > 0, the results of the 10ms 
sweep are already quite close to the long time asymptotic 
state. However, as can be anticipated by the large insta- 
bility growth and asymptotic approach times for w_ < 0, 
the normal state (JJ.) persists. Despite its ultimate insta- 
bility, the 10ms is insufficient for the instability to grow. 
In contrast, with a sweep duration of 200ms, one can see 
the instability of the normal state in this region, although 



= cjQ + C^I'/'P, ip + Kip = —2igr cos 9 . 



(23) 



This pair of coupled first order equations describes the 
exact dynamics of the persistent oscillations. Since the 
motion is in a two-dimensional plane, the attractor is a 
simple limit cycle [t^. In Eq. ([23)) . the phase angle 9 con- 
tinually increases, but has alternate fast and slow regions; 
the motion is faster when \ip\ is larger as may be seen in 
Fig. [TH Such behavior is analogous to a damped driven 
pendulum. In fact, for k ^ luo + U\ip\'^, one may adia- 
batically eliminate \ip\ to obtain 9 ~ {ujq + X) + Acos(20) 
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FIG. 14. (color online). Persistent oscillations at a; = lOMHz, 
{/Af = +40MHz, gyji^ = IMHz starting close to the normal 
state (4) with Sx = Sy = \fN . The top panel shows the at- 
traction towards persistent oscillations, illustrating the tran- 
sient behavior at short times, and the persistent oscillations 
at later times. The inset shows the same data on the Bloch 
sphere using the same shading scheme. For all times after 
~ 15ms, the trajectory on the Bloch sphere is restricted to a 
circle at constant polar angle. The lower panel shows the time 
dependence of l-i/jp and S in the persistent oscillation regime. 



FIG. 13. (color online). Dynamical phase diagram as a func- 
tion of and g ~ g' for U > 0. The phase boundaries at 
Lj± — separate with increasing U and a regime of persistent 
oscillations emerges. The dynamics in this regime is shown in 
Fig. [141 In contrast to Fig. (4] the SRB phase is absent here 
since U > 0. 

where A = U g^r^ /2k. This is the equation of motion for 
a damped driven pendulum, and since > it is driven 
above the threshold required for persistent oscillations. 

B. Phase Diagram for g ^ g' 

Up until now, we have mainly restricted our discussion 
to the experimentally realized case where one necessarily 
has g = g' p^ . However, there are important reasons to 
explore what happens when this condition is relaxed. In 
particular, there are a number of phase boundaries in the 
extended g, g' parameter space, and these can be rather 
close to the experimental situation with g = g' . As high- 
lighted in Ref . [H^l , proximity to these phase boundaries 
is instrumental in explaining regions of slow decay in the 
g = a ' dynamics. Secondly, the proposal of Dimer et 
al (3J| considers a Raman scheme, rather than Rayleigh 
scattering, and involves different hyperfine atomic states. 



In this setting, separate tuning of g and g' could be 
acheived by u sing circularly polarized pump beams in 
a ring cavity |34{ . For these reasons, we consider the 
behavior for g ^ g' . 

In Fig. [15] we set UN = — 40MIIz, and explore defor- 
mations by 5g = g' — g at four different values of fixed 
g = i(g + g'). There are three key aspects to note. The 
first concerns the existence of non-trivial phase bound- 
aries in proximity to the Sg = 0, or g = g' axis. In partic- 
ular, as one transits along the 5g = axis in Fig. I15f d) 
there are two distinct scenarios depending on whether 
\bj\ > l^ttl or < where w„ = UN/2. In the for- 
mer case there is a proximate phase boundary for small 
5g/g, and associated critical slowing down. In the lat- 
ter case the closest phase boundary is horizontal and is 
therefore not crossed by changing 5g. Therefore, for a 
broad range of |a;| < |w„| one may avoid close proximity 
to a phase boundary and the associated critical slow- 
ing down. This w-dcpendcncc of the emergent timescales 
is confirmed in Fig. [71 The second notable feature in 
Fig. [rsl d) is that the SRA and SRB phases, which are 
distinct for g ~ g' , are continuously connected for g ^ g' . 
This may be traced to the lack of factorization of the 
equation of motion for Sj, which simplifies at g = g' to 
S-^ = -ig{ip + ij*){S+ - S-). The third notable feature 
is that there are again regions of persistent oscillations. 
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FIG. 15. (color online). Dynamical phase diagram sXUN = 
— 40MHz, as a function of 5g = g' — g and to for a number of 
values of g — |(<7 + .g')- Vertical dashed lines indicate cuts 
shown in Fig. 1181 



shown in white, where no stable fixed points exist. As 
shown in Fig. 1161 the detailed dynamics differs from the 
persistent oscillations discussed in Sec. IVI Al as may be 
seen from the feature that is no longer constant. 

In general it may be difficult to gain a purely analytic 
handle on the equations of motion in Eq. ^ when g ^ g' 
and U is present. However, in the adiabatic limit with 
K — > oo one may eliminate the photons and consider the 
dynamics of the spins alone. In this limit the equations 
of motion reduce to the following form 
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FIG. 16. (color online). Persistent oscillations for UN = 
-40MHz, ijj = 40MHz, gy/N = 0.998MHz and g'y/N = 
1.002MHz, corresponding to 5g/g = 0.004 or g' /g = 1.004. 
The panels mirror those in Fig. 1141 The shading in the top 
panel match those in the inset, highlighting the initial and 
intermediate trajectories, and the final persistent oscillations. 



where the effective Hamiltonian is given by 

Co{G+Sl + G^Sl) 
H — ojqSz 9,-9 ' 

with u) = uj + USz- The additional contributions 



(25) 



Di = 



Do 



2KrS X (S X z) 



(26) 



(k2 



■S X z, 



S = {S,i/}-Di(S)-D2(S), 



(24) 



are damping terms with G± = {g' ± g)^ and F = g'^ — g^. 
The existence of the photon leakage k therefore means 
that there are two non-Hamiltonian terms in the effec- 
tive spin dynamics. In the special case where U — 
the term Di has the same form as the damping in the 
Landau-Lifshitz-Gilbert equations [t^ [t^ and D2 = 0. 
The former tries to drive the system toward either the 
normal or inverted states. In this U = limit the Hamil- 
tonian contribution in Eq. ([^5]) also reduces to an effec- 
tive Lipkin-Meshkov-Glick (FMG) Hamiltonian frzl - isoj . 

A notable aspect of the adiabatic limit with k — >■ (X), 
is that the resulting dynamics in Eq. (|24)) resides solely 
within the two-dimensional surface of the Bloch sphere. 
The Poincare-Bendixson theorem [t^ therefore excludes 
the possibility of chaotic attractors. In contrast, for 
K ~ 0, the conservative Dicke model with g = g' and 
?7 = is known to exhibit chaotic behavior in the su- 
pcrradiant regime [sil. [3^. In view of this difference, it 
would be interesting to explore the possibility of strange 
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attractors for intermediate k. However, for the param- 
eters we have explored numericaUy we see no evidence 
for strange attractors. Indeed for g — g' and U > 
we have demonstrated the existence of a hmit cycle gov- 
erned by Eq. ((23| . Nonetheless, it is worth noting that 
the nonlinear equations of motion in Eq. ([3]) are cl osely 
related to the Maxwell-Bloch equations for a laser |8l| . 
These are known to be equivalent to the paradigmatic 
Lorenz equations (82j , the archetypal example of dissipa- 
tive chaos. However, an important difference from the 
Maxwell-Bloch equations is the absence of external driv- 
ing in Eq. ([3]). It would be instructive to explore the 
ramifications of this in more detail. For work on chaos 
in a closely related optomechanical system see Ref. [s^ . 
Further discussion of the phase diagram for g ^ g' is 
given in Appendix lEl 



VII. BEYOND THE EFFECTIVE DICKE 
MODEL AND ITS SEMICLASSICAL 
TREATMENT 

In the preceding sections, we have discussed the semi- 
classical dynamics of the nonequilibrium Dicke model, 
and its relation to experiments on the self-organization 
of BECs in optical cavities [l^ . Within the semiclassical 
description of this model we have found a rich variety 
of stable attractors including non-trivial steady states, 
persistent oscillations and regimes of bistability. Having 
established a wide variety of predictions for the semiclas- 
sical behavior of the open Dicke model, we now consider 
what effects may arise in going beyond this effective de- 
scription. In Sec. IVII Al we first consider modifications 
to the Dicke model itself, arising from higher momentum 
states and other terms in the effective Hamiltonian. In 
Sec. IVII B] we briefly comment on the possible modifica- 
tions due to higher order quantum effects. 

A. Modifications of tlie Effective Dicke Model 

As outlined in Sec. [Hi the derivation of the effective 
Dicke model involves a projection onto the subspace of 
the two lowest lying momentum states; see Appendix VK\ 
Without this projection, there would also exist coupling 
to higher momentum states such as 




In general, the occupation of these excited states is ex- 
pected to be small for low intensity cavity light fields, 
as supported by the time of fiight images of Ref. p^ . 
Nonetheless, there are regimes of parameter space where 
these states may be important. In particular, these high 
momentum states may destabilize certain phases pre- 
dicted by the reduced Dicke model. Specifically, the 
inverted state involves excitation to the north pole of 



the Bloch sphere and may be susceptible to destabiliza- 
tion. Indeed, in the parameter regimes where the ef- 
fective Dicke model predicts the inverted state, kinetic 
approaches predict heating [s^] ■ Likewise, regions of the 
superradiant phase in which the majority of the atoms 
are in the non-zero momentum state may be unstable. 

Although the stability of these particular states may 
be modified, we anticipate that many of our predictions 
are only weakly affected by these additional states. This 
is supported by the clear quantitative agreement between 
the experimentally observed onset of superradiance and 
the reduced Dicke model [22|. Our findings within the 
projected subspace may also describe experiments ex- 
ploiting internal hyperfine states (33 |. where no higher 
levels exist. 

In addition to the effects of higher momentum states, 
one should also note that the intensities of the forward 
and backward propagating pump beams are not of equal 
magnitude in the experiments of Ref. ; they differ by a 
factor of 0.6 due to losses on reflection. This introduces a 
coupling to the state ^ ^^_|_ /3|afc, j3k) which exhibits 
odd parity under reflection in the pump direction. We 
expect such contributions to play a similar role to higher 
momentum states. 

A further source of possible departure from the ide- 
alized Dicke model with g ~ g' arises because of the 
finite atomic recoil energy. As we demonstrate in Ap- 
pendix [XJ] the processes leading to g and g' correspond 
to different detunings of the intermediate states, so that 
igjg = Wr('j-'c — a;p)/2(a;a — Wp). We have investigated the 
possible impact oi g ^ g' in Sec. IVI Bl and in Appendix 
IeI where we showed that differences 5g/g of the order of 
10"'^ can cause one to cross a phase boundary. However, 
for the typical values of the detunings used in Ref. (2^ . 
this asymmetry is of order 5g/g ^ 10~^^. It is therefore 
too small to have any significant effect on our findings. 
Nonetheless, in experiments with a smaller atom-pump 
detuning, this asymmetry could play a crucial role. 



B. Corrections to Semiclassical Dynamics 

In the above discussion we have focused primarily on 
the semiclassical dynamics of the effective Dicke model 
through the solutions of the equations of motion given 
in Eq. . In some places we have also incorporated the 
leading effect of quantum fiuctuations by using Wigner 
distributed initial conditions. This may be interpreted 
as including subleading 1/A'^ corrections where N is the 
number of atoms [69j . More generally it would be prof- 
itable to investigate the full quantum dynamics governed 
by the density matrix equations of motion in Eq. ([T|). 
However, it is important to bear in mind that since the 
density matrix describes an ensemble average it will in 
general mask the effects of spontaneous symmetry break- 
ing. It will also wash out collective and persistent oscilla- 
tions (ssf . Nonetheless, in a single experimental run one 
still observes spontaneous symmetry breaking (22l . [23j . 
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and the density matrix describes the average over many 
runs; see for example Ref. [s^. In order to recover in- 
formation on these non-trivial features within the density 
matrix formulation, one may consider higher order corre- 
lation functions. We leave this problem for future work. 



Appendix A: Derivation of the Generalized Dicke 
Hamiltonian 



Here we provide a derivation of the generalized Dicke 
Hamiltonian for ultra cold atoms placed in an optical cav- 
ity, as illustrated in Fig. [T] For simplicity, we consider a 
homogeneous system, and neglect the effects of the finite 
beam waist of the pump and cavity fields: 



VIII. CONCLUSIONS 

In this manuscript we have explored the collective dy- 
namics of ultra cold atoms in transversely pumped op- 
tical cavities. Within the framework of the effective 
noncquilibrium Dicke model we present a detailed discus- 
sion of the rich phase diagram of asymptotic attractors, 
including steady states, coexistence phases and regimes 
of persistent oscillations. We show that the inherent 
timescales for the destablization of the initial state, and 
the decay time towards the asymptotic attractors, show 
strong variations throughout the dynamical phase dia- 
gram. Crucially, we have demonstrated that two distinct 
principal timescales emerge, corresponding to the energy 
scales ojo and ujq/k- The scale ujq, characterizes both the 
typical frequency of collective oscillations and their decay 
rate for a broad range of parameters. The slower scale 
Wq/k, governs the decay rate in proximity to dynami- 
cal phase boundaries, and may be interpreted as critical 
slowing down. Most notably, in the regime a; < w„ , sweep 
experiments over 200ms may be required in order to reach 
the asymptotic regime. It would be profitable to explore 
this experimentally and we discuss the broad implica- 
tions for finite duration experiments. In particular the 
superradiant phase divides into two distinct regimes, de- 
noted SRA and SRB, with the relaxation rates ujq/k and 
Wo respectively. From a theoretical perspective it would 
be valuable to investigate the role of quantum fiuctua- 
tions, and the effects of states outside of the two-level 
Dicke model description. It would also be interesting to 
explore the ramifications of these findings in other real- 
izations of the Dicke model. 
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goitp + ip^)cos{kxi) 



-1-0/ cos{kzi — cupt) + fib cos{kzi + cupt) 



(Al) 



The Hamiltonian acts on both the centre of mass posi- 
tion of the atoms, and their electronic state. The lat- 
ter is restricted to the two states involved in the opti- 
cal transitions, denoted by e and g for the excited and 
ground states respectively. In this basis a^'^ = |e)i(e|i 
and (T,^^ = |c)i(<7|i. The sum over i runs over the number 
of atoms present in the cavity, N. The matter-light inter- 
actions correspond to the dipole coupling of the atomic 
transition to the fields of the cavity and the forward and 
backward pump fields. The cavity field ip has explicit 
quantum dynamics, while the time dependence of the 
pump fields is externally imposed. The cavity-atom cou- 
pling is designated go. The pump strengths fif{b) de- 
scribe the pump beam in the forward (backward) direc- 
tion, where we allow for imperfect retro-reflection as dis- 
cussed in Ref. [l^] and Sec. I VIII We have neglected any 
difference between pump and cavity wavevectors. 

The matter-light coupling in Eq. (|A1[) contains both 
CO- and counter-rotating terms. The rotating wave ap- 
proximation consists of neglecting the counter-rotating 
terms on the basis that the detuning is large compared 
to the coupling strengths. This approximation is valid 
here, since Wc, Wq, Wp are all optical frequencies, and of 
order 400THz. This is to be compared to the coupling 
strengths flf^bV^ goVN ^ 1 GHz. Working in the 
rotating frame at the pump frequency ujp, and neglecting 
the counter-rotating terms, the Hamiltonian (|Aip can be 
rewritten in the form H = Hq + Hi where 



Ho = Ac^^V + Aa 51 



— y v2 



(A2) 



and 



Hi = 



Of -, 

go^' cos.{kxi) ^ -l-e"-^- 



O 



H.C.. 



(A3) 

Here Ac = Wc — is the cavity-pump detuning and 
Aq = cJq — cjp is the atom-pump detuning. 

As described in Sec. [Hi the effective Dicke model is 
a low energy description within the electronic ground 
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state, valid if the coupling to the electronic excited state 
is small compared to the detuning Aa- As such, one may 
proceed by making a Schrieffer- Wolff [s^l transformation 
and eliminating the excited electronic state. This gives 
a transformed Hamiltonian H ~ + {i/2)[S, Hi] where 
[S, Ho] = iHi . One should choose S so that H has no 
coupling between the resulting dressed electronic states: 



5 X 10 In this approximation Eq 
1 



iS 



where for simplicity, position and momenta are now ex- 
pressed in units of the cavity wavelength. In addition, 
wc have suppressed the atom labels and the summation. 
Here the differential operators / and g are given by 



1 



where uir = h^k^ /2m is the recoil energy. The resulting 
Hamiltonian has the form 



H 



4A„ 

„2 



4]) becomes 

[nj + nl + 2nfnb cos(2z)] 



A, 



■ cos^(a;) 



go 
4 



1 

a: 



cr'^f-H.C. 



A. -A '""'^^^ 
iiflf - rJb)sin(z)(V'^ - V') 



(A5) 



We proceed by introducing a basis set of atomic center 
of mass states for the atoms in their electronic ground 
state. The first two of these states will correspond to the 
"spin down" and "spin up" states of the effective Dicke 
model. These basis states are given by the eigenstates of 
the first line of Eq. (|A5P , which may be written as 



f ' (A6) 



The energies are given by E' = a;r(?^+ecr,m)+const. where 



H = Ho-^f~^r 



r2 f^i) 



- f V'V [5(1 + + ,r (1 + e^n + h.c] 



lAH^e'^ + ^be-")ige-''' + .g*e") + H.c. 



^)(%/e-^^+06re'-)+H.c., (A4) 



where /* and g* denote the complex conjugates of the 
Hermitian operators / and g. In writing Eq. (jA4| we 
have eliminated the excited electronic states, but have 
made no further approximations. As a result, the form 
of H is rather unwieldy. In order to expose the resulting 
behavior wc will consider two classes of approximation: 
the small recoil approximation ujr/ Aa ^ 1 and the weak 
pump approximation q = £l/f2f,/(4a;r.Aa) <C 1. Both 
of these approximations are valid for the experiments of 
Ref. [23| and we will now consider each in turn. 



1. Small Recoil Approximation LOr Aa 



If the recoil energy a;^ is small compared to the atom- 
pump detuning Aa , then one may set cj^ = as a first ap- 
proximation. In this case / ~ A"""^ and g ~ {Aa — Ac)^^ 
become c-numbers and the form of Eq. (|A4p simplifies 
considerably. This ap prox imation is well justified for 
the parameters of Ref. [23 as w^/Aq ~ 50kHz/lTHz ~ 



+ So 



2gcos(2z)] (j)a 







(A7) 



is the Mathieu equation [88| and cr = ± label the 
even and odd solutions. The Mathieu parameter q = 
n fill, / (AAaUJr) is a dimensionless measure of the pump- 
ing strength. Due to the form of the matrix elements aris- 
ing from the terms linearly dependent on g^ in Eq. (|A5p , 
not all of the configurations of cr, m, n are coupled. Only 
those states that can be reached by a sequence of absorp- 
tion and emission processes, starting from an atom in the 
ground state, need to be included. If fi;, = 51/, then only 
the even Mathieu functions </>+.,„ (z) need to be included. 
However, if 51;, 7^ 51/, the odd Mathieu functions (/)_., „(z) 
should also be considered. 

To recover the effective Dicke model one must restrict 
attention to the two lowest states, and work in the limit 
ilf ~ rii, — rt. In this case, the lowest coupled states are 



and $ 



+,1,1 



The values of the parameters oj,ujo:U 



in the effective Dicke model in Eq. 1^ can be found by 
evaluating ($+,o,o|-ff |$+,o,o) and In 
terms of the Dicke model parameters in Eq. ^ , the ener- 
gies of a configuration with riph photons and all N atoms 
in either their ground or excited states are given by: 



UN\ 



(A8) 



By comparing and with the expressions for 

(<I>+,o.o|^|*+,o,o> and (<I>+,i,i|i/|$+4,i) one may iden- 
tify the coefficients uj, loq and U. We find 



Ar 



5glN 



8(A, - A,) 



U 



9l 



^{Aa-A,y 



(A9) 
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where we have carried out the summation over atoms and 
made use of the results ($+.0,0 1 cos^(a;)|$+.o,o) =1/2 and 
cos^(a;)|<i>-|-,ij_) — 3/4. These cocfEcients agree 
with those of Ref. |23l when the pump and cavity fre- 
quencies are near detuned. The two-level energy splitting 
is given by the difference of the eigenvalues of the states 
written in Eq. (|A6p and so 



.0). 



(AlO) 



Evaluating the off-diagonal elements ($+_o.o|^|*i'+.i,i) 

and equating ($+_o,o|^|*&+a,i) ~ ffV'^ + .9^ one finds 
the remaining Dicke model parameters 



9 = 9 



90^ 



1 

a: 



1 



dz (f>^_Q{z) cos{z)(f)^_i{z). (All) 



Up until this point all the results we have derived in 
Eqs. (|A9[) - (|A11[) are formally exact for arbitrary q = 
ri/r2b/(4AoWr)- However, the fact that one can restrict 
to the two lowest momentum states is only valid for weak 
pumping, i.e. small q. We will therefore focus on the 
small q limit; in Sec. lA 31 we will return to the general 
q case, including also the presence of higher momentum 
states. In the small q limit one obtains 



UJQ « 2uJr, 9 = 9' 



90^ / I 1 
4 Ia„ A„-A, 



, (A12) 



where we have used the approximations (j)^^a{z) « l/v27r 
and sa cos(z)/-y/7r. If one further neglects the cavity- 
pump detuning Ac in comparison to the atom-pump de- 
tuning Ag, then these expressions reduce to those given 
in Ref. [23 • In Sec. IA2l we will generalize the results of 
this section to include the effects of nonzero ujr/Aa- 



2. Corrections Due to Non-Zero uur/Aa 

In order to quantify the effects of non-zero w^/Aa, we 
return to Eq. (|A4p . but continue to make use of the q = 
approximation used in the second half of the previous 
section. The main difference that non-zero LUr introduces 
is that g and g' are no longer equal. One obtains 



9 
9' 



90^ 
4 

90^ 
4 



1 



1 



Aa + UJr - Ac Aa - LLJr 
1 1 



Aa - LOr - Ac Ag + Ur 



(A13) 



where for simplicity we set f2/ = 51b = SI. In the limit 
-J' 0, these reduce to Eq. (|A11|) . In Sec. IVIBI we 
show that a phase boundary can be crossed if the frac- 
tional difference 5g/g is large enough. At leading order 
in Wr/Ao the fractional difference given by Eq. (|A13[) 
is 6g/g = w^Ac/A^. For the experimental parameters in 



Ref. [221 this fractional difference is too small to cross the 
phase boundary. However, for smaller A^ this fractional 
difference may become significant. 

The results in Eq. (|A13p are obtained by the 
same procedure as in the previous section, by eval- 
uating the off-diagonal matrix element and equating 
($+,o,o|-ff|'&+4,i) =5V'^+5'V'- In deriving Eq. (|A13|) we 
use the basis states ^+fifi{x, z) = l/(27r), 2;) = 

cos(x) cos(2:)/7r, which follow from Eq. (IA7|) at g = 0. We 
further employ the identities 



dC 



cos(C) = -, 

X 



cos(C) 



X — IbJr 



where V ~\x — i2ijjy.dc\ ^ . 

The remaining parameters of the Dicke model are 



9l 



1 



(5 + 3a;,. 5 — uj,. 5 + uJr 
1 2 



9l 



5 + 3a;,. 5 — UJ,. 5 -\- u:r 
1 2 



Aa + 3a;r A^ — a;^ Aa + w,. 



where 5 = Aa — Ac = a;^ — Wc. In the limit of a;,- — )■ 
these reduce to Eq. (|A9p and the small q expansion 
of Eq. (jAlOp . These expressions are found using the 
same procedure as outlined in the previous section, by 
equating E^,E^ in Eq. (jA8|) with the expressions for 

($+,o,ol^l$+,o,o) and To evaluate 

these expressions we use the results /|<E'-)-.o,o) = (Aq + 
a;,.)-i|$+,o,o) and g|4>+,o,o) = (A^ + a;,. - Ac^^ |$+,o,o) 
together with 



dC 



cos(C) 



(l + e-2<)+c.c.+ 



x — i2ujrd(; 
c.c. 



(1 



„2< 



X — i2LJrd( 



where x = Aa + ujr and x = Aa 
involving / and g respectively. 



X — 2a;.r X -|- 2a;r 
— Ac for integrals 



3. Equations of motion in extended state space 

When the parameter q is not small, the restriction to 
two atomic states is no longer valid. Including higher 
momentum states, it is no longer possible to map the 
problem on to an effective spin Hamiltonian. However, 
it is still possible to derive a semiclassical description 
of the coupled atom and cavity system. In this gen- 
eralised case, the semiclassical description consists of a 
Gross-Pitaevskii equation for a macroscopically occupied 
atomic wavefunction x(a;, z) coupled to a quasi-classical 
Heiscnberg equation for the photon field V'- 
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Decomposing the atomic wavefunction in some basis 
xix,z) = J2aXa^aix, z) and splitting the Hamihonian 
into the parts 

H = + V-^/i^^^ + (V^t + ^,)/i(2) + j(^t _ ^)/j(3) 
the exphcit equations of motion read 



(A14) 



(A15) 



where have defined Mj^^ = ($q,|/i^"-'|$^). As long as q is 
small, one may truncate to two atomic basis states, and 
thereby recover the semiclassical equations in Eq. ([3]), 



where S, 



IXil 



Ixol and + iSy 



XiXo- 



Appendix B: Fixed Points with Arbitrary g and g' 

In general it is difficult to obtain explicit closed form 
expressions for the steady state solutions of Eq. ([3]). How- 
ever, in the special case where U = simplifications occur 
for arbitrary g and g' . More generally, for arbitrary U, g 
and g' one may obtain self-consistent implicit solutions. 
We discuss these cases below. 



1. f/ = 

In the case where U = the nonlinear equations be- 
come linear equations for the variables tpi, -02, Sx and 
Sy, where Sz enters via the coefficients. Decomposing 
Eqs. © and ^ into their real and imaginary parts yields 



-ujoSy = 2{g - g')Szi>2, 
KTpi - uj-tjj2 = -(.g - g')Sy, 
ujipi + K-tJj2 = -(.g + g')Sj;. 



(Bl) 



The last two equations may also be written in the form 

(w^ -f- K'^)i>i = -u){g + g')Sx - K(.g - g')Sy, 



(w^ + K^)'ip2 = -K{g + g')Sx + uj{g - g')Sy. 



(B2) 



The condition for non-trivial solutions yields the dctcr- 
minantal self-consistency equation 

^{g^~g'^fSl+AujLOo{g^+g'^)SM^^+^^)^l^^- (B3) 
This may be solved to yield 

-ujuonig^ + g'^) ± ^ {2uL0^gg'Y - Loln^ig^ - g'^f 



S 



corresponding to a non-trivial superradiant phase with 
i/' 7^ 0. In the limit g = g' one obtains Sz = —loq{uP' + 
K^)/8ajg-^. The critical coupling strength for the onset of 
superradiance corresponds to Sz = —N/2 and is given by 



g^/N 



4a; 



(B5) 



in agreement with the results of Dimer et al [3J|. It 
also coincides with Eq. (|lip for the onset of the SRA 
phase when [/ = 0. The explicit dependence on k in 
Eq. (jBSP emphasizes that the transition occurs in an open 
system. In the limit k = one recovers the location of 
the superradiance transition, g\/N = ^^lijujq/2, for the 
equilibrium Dicke model with counter-rotating terms. 



2. [/ / 

In the general case with arbitrary U, g and g' the 
steady state solutions of Eq. ^ are more difficult to ob- 
tain in an explicit form. However, one may still obtain 
self-consistent solutions which relate the photon density 
to Sz for example. In turn, these implicit consistency 
equations may be solved numerically. For generic pa- 
rameters the steady state equations of motion may be ob- 
tained from Eqs. (jBip by replacing loq — >■ uq and — ^ a), 
where u) = lu+USz, i^o = loq+Uu and n = ~ ■ipi+ip2- 
For a given photon occupation number, n, non-trivial so- 
lutions satisfy the determinantal Eq. (jB3[) with these re- 
placements. Using the explicit form oi uj = lj + USz one 
obtains a modified quadratic equation for Sz'- 

(x' - l6g^9'^)S^_ + 2loCjoxSz + [uj^ + = 0, (B6) 

where we define x = 2((7'^-|-(7'^) + ?7a)o- This has solutions 



Sz 



Wo 



x^ - ^Qg^g 



(B7) 

In order to have real solutions to Eq. (jB7[) one requires 
16g^g' (w^ + K^) — K^x^ > 0. This translates into the 
condition that max(0, n_) < n < n-|_ where 



?i+ = u- 



-2(.g2 + g'^) - Uujo ± igg'n-WuJ^ + 



One should restrict attention to those cases where l&j < 
A^/2. Having found a solution for Sz in terms of the 
number of photons n, one may also find an equation for 
n in terms of Sz- Using the analogues of the first two 
equations in Eq. (|B1[) with luq — !■ ljq 



Si 



1 



Hg + g'r ^g-g'y 



(B8) 

where r = Sy/Sx- Using the fixed length spin constraint 
g2 ^ ^2/4 Qjjg j^j^y eliminate S^ = {N^ /i- S^)/il + r^) 
(B4) in favour of the ratio r. This ratio may be obtained by 
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using the analogues of the first two of Eqs. (jBl[) to sub- 
situte tpi and '02 into the analogue of the third equation: 



This yields a quartic equation 



9 + 9 ' 



LdQK, 



Substituting into Eq. (jB8|) yields 



(B9) 



(BIO) 



Appendix C: Linear Stability for Arbitrary States 



In Sec. IIII CI we linearized the equations of motion 
around the normal (^) and inverted (ff) states. The cor- 
responding self-consistency equation is given by 



0. (CI) 
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/ Cj — IK 



-2gS^^ + Ur^S^ 
2g'S^o - UroS+ 
\ 9So-g'So 





— (cj + zk) 
-2g'S^ + Ui^oSo 
2gS^-U7jjoS^ 
-9S0 +9'S^ 



= [(77 
T277(77 -f 



-f Zk)^ - L 

iK){g^N 



]) + {g''N-g''Nr 



g'^N)T2{g^N + g'^N)u;^u;o, 



(C2) 

whose roots characterize the possible instabilities. In the 
case where g ~ g' Eq. (|C2p reduces to Eq. ([T5|) . More 
generally one may linearize the equations of motion in 
Eq. dSl) around an arbitrary state, so that ip = + Sip, 
= Sq + 6S- and S'' ^ + SS'', one obtains 

5S^ = — iCjqSS^ — ill {iPqStP + 4'aSip*)SQ 

+ 2z(57^o + g'ro)SS' + 2t{g5i, + g'Sr)S'o. 

SS' = - igSi'S+ + igSi'*So + ig'SipS^ - ig'5i'*S+ 

— igtpoSS^ + igt/jQSS^ + ig't/joSS^ — ig' iI'qS S'^ , 
Stp ^ — (k + iu})S'ip — iUipoSS^ — igSS^ — ig'SS'^ , 

(C3) 

where uj = oj + USq and ujq = ujq + U\tpo\^. Parameter- 
izing 6^ = ae-"'* + 6*e"'**, dS' = ce"*''* + d*e''''* and 
SS^ = /e"*''* -I- /*e"'**, one obtains a set of algebraic 
equations for the coefhcients a, b, c, d and /. The corre- 



sponding secular equation is given by \rjl — ] 
I is a 5 X 5 unit matrix and 



0, where 



9 9 

-9' -9 

cjo 

-Wo 

-(.gV'o + .^Vo) (5^0 + g'%) 



-uro 

-2(300 + 5VS) 
2(50* +gVo) 




(C4) 



In general, there are in fact only four independent equa- 
tions owing to the fixed length constraint, = N'^/A. 
This is reflected by a redundant zero mode, 77 = 0, cor- 
responding to longitudinal fluctuations in the length of 
the spin. Although Eq. (jC4|) captures all of the essen- 
tial information, it is convenient to eliminate this zero 
mode and use a 4 x 4 matrix representation for the 
physical transverse degrees of freedom. Differentiating 



the fixed length constraint with respect to time yields 
25^5*^ + [S+S- + S-S+) = 0. Linearizing around a 
fixed point gives 2S^6S^ + S+SS~ + SqSS+ = 0, as 
may verified by Eq. (jC3p and Eq. ([3]). Using the nor- 
mal mode parameterization one obtains the relationship 
/ = —{S^c + S^d)/2SQ, between the coefficients. Elim- 
inating / from the linear equations yields |7;I — M| = 0, 
where I is a 4 x 4 unit matrix and 



UJ — IK 



-2gS^^ + Ui'lS^ 
V 2g'S^^ - (700* 5o+ 





— (w + in) 
-2g'S^ + C/0o5o- 
2,95o^ - [/7Ao5o+ 



Wo 



g - Utlj„S+/{2S^) 
-g' + Ur„S+/{2S^) 
(g0o+ffVo*)'5o+/5o^ 



[/0o5o-/(255-) 
-C/0o*5o"/(2^o) 
(5^o+5Vo*)5o-/^o' 



(C5) 



-(300* + 9'i'o)S+IS^^ -wo - (500 + g'i^o)S^/S^J 



1. Stability of the normal (ij.) and inverted {■fi) 
states with arbitrary g and g' 

In the case of the normal (^) and inverted (ff) states 
where 0o = and = ^N/2, the determinantal equa- 



tion |7/I — M| =0 reduces to Eq. (|C1|) . This corresponds 
to the quartic equation given in Eq. ()C2p . In order to 
find when the roots become unstable it is convenient to 
decompose the quartic polynomial into its real and imag- 
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inary parts so that Azf: (77) + iB^ {rj) = where 

= {rf - 4 - K'){rj' c^) + (g'-N g'^ Nf 

(C6) 

We may thus find the roots of the equation B^{r]) = 0, 
and subsequently impose the condition A-^{if) = on 
these solutions. It is readily seen that Bz^ijj) = has a 
solution 77 = corresponding to an exponentially growing 
mode without oscillations. Substituting this value into 
the expression for A^{rj) yields the condition A^{Q) — 0: 

(4 + ^2)^2 + (g2^ _ g,2^^2 ^ 2(g2 TV + g'^N)u:^u:o - 0. 

,(C7) 

In the special case where g — g' one recovers the critical 
condition for the onset of the SRA phase as given by 
Eq. pT|) . Alternatively, there is also an instability with a 
finite frequency 77^ 



Jq ± (f/^iV — g' N) corresponding to 
B^{rj) = 0. Demanding that A^{if) = yields a critical 
condition for the ratio of the couplings 



q' _ (w=p + ujq) + K 



2 ' 



(C8) 



rather than their absolute scales. This instability condi- 
tion manifests itself as the phase boundaries shown in the 
upper panels of Figs. [12] and [TOl Denoting Sg = g' — g, 
g = {g+g')/2, and noting that ujq ^ n, the phase bound- 
ary given by Eq. (jC8[) may be approximated as 



(C9) 



dg 2ajoWzp 
g oj'^ + k'^' 

In a similar fashion Eq. (jC7p may be recast as 

1 — Sg^N/idoUzfi ' 



(CIO) 



where g J is given by Eq. (|lip . For the range of parame- 
ters shown in Figs. [T5] and [TOl Sg < g,gf, and Eq. (|C10|) 
is effectively independent of 5g/g; this yields the horizon- 
tal phase boundaries in Figs. fT5l and [T9l 
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FIG. 17. (color online). Magnified portion of the bottom 
panel of Fig. [4] in the vicinity of the tricritical point where 
three phase boundaries cross. In addition to the phases visible 
in Fig.|4]there is a narrow region denoted as 2SRA, where the 
two distinct SRA solutions given by Eq. (|10p coexist; see inset. 



respectively. The eigenvalues satisfy Eq. (PT|) which is 
the square of a quadratic equation. The exact eigenval- 
ues corresponding to fluctuations around the stable SRB 
fixed points are given by Eq. ([22]) . 



Appendix D: Transitions Near the Tricritical Points 

As noted in Rcf. [5(| for UN < -2k, three of the 
phase boundaries in Fig. |4] cross at the point w = 
y^oj^ — K^, g = a;o?7/4. As shown in Fig. [iTl in this 
vicinity there is a narrow region where the two distinct 
SRA solutions given by Eq. pUj) . together with their par- 
ity symmetry partners, are stabilized; sec also Fig. [5] On 
Fig. m these occur within the width of the line marking 
the boundary of the SRA -I- -ff and SRB regions. At g = gt 
the two pairs of SRA solutions merge, and switch to two 
pairs of SRB solutions. After this, one of each pair is sta- 
ble whilst the others are unstable, as generically occurs 
in the SRB phase. 



2. Stability of SRA and SRB with g = g' 



Appendix E: Further cuts through the phase 
diagram with g ^ g' and [/ 7^ 



When 5 = 5' the SRA phase has 5^ = and 5*^ = 5^. 
In this case the matrix M has eigenvalues rj, satisfying 



''-^\2gS^,-U^^,S^,\^ 

Jr. 



0, 



(Cll) 



where uj = oj + U Sq and cjq = + C/|V'oP- 

When g = g' the SRB phase has uj ^ ujq ^ and 
= ipi + i'ijj2 is purely imaginary. In this case the di- 
agonal blocks of M are proportional to unity and zero 



In Sec. IVI Bl we presented the dynamical phase dia- 
gram with g ^ g' and UN = — 40MIIz, illustrating the 
dynamical phase boundaries which emerged for small dif- 
ferences between g and g' . Here we provide further cuts 
through the phase diagram in order to fully expose the 
rich topology. In Fig.[TS]we show a sequence of cuts with 
UN = — 40MHz for different values of g' /g. These may 
be compared with the bottom panel of Fig. (|4]) which 
has g' /g = 1. In view of the duality relation in Eq. ([5]) 
we only show the results for g > g' . For complete- 
ness, in Figs. [12] and [501 we also show cuts of constant 
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FIG. 18. (color online). Evolution of the phase diagram 
shown in the bottom panel of Fig. 2] as one varies the ra- 
tio g' /g, increasing from 1.000 (top) to 1.004 (bottom). We 
use the same phase labelling conventions as in Fig. 21 Vertical 
dashed lines indicate cuts shown in Fig. 1151 



g = + g')/2 and g' / g respectively, with U = +40MHz. 
The central white region in these figures is continuously 
connected to the regime of persistent oscillations de- 
scribed in Sec. IVI Al We note that the regions where 
the normal and inverted states are stable have identi- 
cal shapes to those seen for U = — 40MHz in Fig. [T5] and 
in Fig. [THl respectively, but are displaced vertically. 
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(c) g\'N=0.9MHz 
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FIG. 19. (color online). Analogue of the phase diagram in 
Fig. [15] for UN = -|-40MHz. We use the same phase labelling 
conventions as in Fig. ID The panels show the dependence 
on Sg = g' — g for fixed g = |(<? -f g')- The white region 
corresponds to a regime of persistent oscillations. Vertical 
dashed lines indicate cuts shown in Fig. (20] 



Appendix F: Wigner Distribution 

Sampling initial conditions from a Wigner distribution 
is readily achieved by combining a Holstcin-PrimakofF 
representation for the collective spin operators, Sz = 
-N/2 + a^a and S' ~ ^Na [sl [90] with a har- 
monic oscillator decomposition of the auxilliary bosons, 
a = (q + ip) I The corresponding Wigner distribution 
Wlq^p) = e~'i ~P /tt reflects the Gaussian ground state 
wavefunction of the Harmonic oscillator. In terms of the 
Bloch sphere coordinates with 5~ ~ sin(6')e~"^ 
this corresponds to ~ A/(g^+p^)/(iV/2) and = 
— Arg((3' + ip). In a similar fashion, the Wigner distribu- 
tion of the initial photon field ip = {Q + iP)/V^ is given 
by W{Q, P) = e~'^ ~^ /tt. In order to employ these en- 
sembles we sample initial conditions from these distribu- 
tions of [q^p) and (Q, P) and time evolve the scmiclassi- 
cal equations of motion. We then average the final results 
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